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ABSTRACT 



This thesis develops a program which will merge or 
overlay imagery and terrain elevation data and create a 
synthetic 3-D perspective view of the ocean bottom. The 
observer may position himself at various locations and see 
the terrain from different viewpoints. The elevation data 
is grouped into triangular panels and the color 
information is averaged from the imagery data file. The 
entire panel is assigned a single color equal to the 
average. These panels are then projected onto an image 
plane by using a 3D to 2D perspective transformation. 
Hidden surfaces are removed by a "painters" algorithm 
which relies on sorting the panels based on distance from 
the observer. 
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INTRODUCTION 



The goal of this thesis is to develop the software 
required to display a 3-D perspective view of combined 
sonar imagery and bathymetric data. The result is a 
synthesized image representing a 3-D perspective view of 
the terrain from a given observer location. 

The issues involved include: first, how can the sonar 
imagery data and the bathymetric data be combined; second, 
what is the methodology of creating the 3-D perspective 
view from all observer locations; third, how fast can each 
view be generated using real data and fourth, how is the 
picture quality affected by the resolution of the data. 

Conventional methods of displaying elevation data are 
with contour lines or 3-D grid line drawings. Some 
applications merge the contour lines with color data and 
this aids in the overall comprehension of the data. 
Recently the speed of the digital computer has been called 
upon to generate 3-D views based on elevation with 
shading. The shading value may be based on the elevation, 
or may be based on other information available. These 
displays are an improvement over the simple contour plots, 
as well as the 3-D grid line drawings. The advantage of 
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any 3-D perspective display method is the ability to 
"see" the terrain in a fashion which appears normal to 
the observer. 

Applications of this method range from low cost flight 
simulators for manned or unmanned vehicles to real time 
displays of multiple source information in a fashion which 
results in a greater level of comprehension or 
interpretation. For example, a combat aircraft pilot must 
be able to assimilate vast amounts of data. The pilot must 
search multiple electronic displays in order to put 
together a picture of the environment in his mind. The 
pilot must answer several questions simultaneously: Where 
am 17 , Where are my friends?, Where is the enemy?. What 
is the condition of my aircraft? Simplifying the 
presentation of this information in a simulated 
perspective view of the surrounding environment will 
improve the response of the pilot. [Ref . 1, p. 64] The 

ability to combine the various data types into a combined 
display is a form of "sensor fusion". In the aircraft 
example the data may be video from external video cameras 
and radar or infrared sensor returns. 

Generally, sensor fusion refers to the ability to 
merge layers of information from various sources into a 
new form. The data may be imagery data from Landsat or 
aircraft, or any' other layer of data within the same 
geographical area such as magnetic anomaly or infrared 
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response data. These various types of data must be in a 
form which can be represented by a specific color or gray 
shade. Adding color or shading from the other data source 
allows the correlation between it and the terrain 
elevation data to become apparent. A specific area of 
interest includes shipboard systems where numerous sensors 
are gathering data and the operator selectively views the 
output from a single sensor or a combination of sensors. 
If the outputs from multiple sensors are combined into a 
3-D perspective view of the surrounding environment, the 
response time of the operator will improve. By utilizing 
the 3-D techniques, a view forward of the ship could be 
displayed on a CRT screen. Combining radar elevation data 
with an image file database would permit the formation of 
a 3-D perspective view of the radar image. This same 
technique used with the sonar data obtained by the ship's 
sonar systems and imagery data files of the harbor channel 
would allow a realistic view of the channel and permit 
more efficient usage of the available sonar systems. 

The software developed in this thesis overlay side- 
looking sonar imagery data of the ocean bottom onto the 
bathymetric data of the same geographic area. This is an 
extension of work by L. Coleman [Ref. 2]. Coleman's work 
concentrated in generating a 3-D perspective view of land 
terrain. Several items added in this work include the 
ability to approach the area from any heading (observer 
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locations within the area boundaries are permitted) , 
reduction of the time to create one image, and drawing of 
the image is included in the main routine. In addition to 
the work by Coleman, a separate project by McGhee, Zyda, 
Smith and Streyle incorporates many of the same techniques 
[Ref. 3]. The imagery data is in the form of a 512 x 512 
pixel image with 256 possible gray levels. The bathymetric 
data is gridded and consists of depth values on a latitude 
longitude grid. The bathymetric data, called the elevation 
data, is segmented into triangular panels. The boundaries 
of each panel overlay the image data, image pixels inside 
the panel are averaged and this average gray level is 
assigned to the entire panel. To generate the perspective 
view, each elevation point is projected onto an image 
plane located at the observer location, and oriented with 
respect to the observer course. This projection is 
accomplished through a 3D~2D perspective projection 
transformation which maps the elevation point coordinates 
to image plane coordinates. The synthesized image is 
generated on a monitor and approximates the view from the 
observer position with a 38.6 degree field of view. The 
synthesized view is directly affected by the resolution of 
the input data. The elevation file resolution affects both 
the elevation drawing and the shading of the image. 
Shading is affected due to the averaging of pixel shades 
within each triangle. The drawing is affected by the low 
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sampling rate of the terrain; only the larger features 
will be present in the displayed image. 

The software is implemented in FORTRAN on a DEC Micro- 
Vax GPX II. This system is capable of 1024 x 864 pixel 
resolution, with 256 colors or 256 gray levels. A graphics 
software package from DEC, MicroVMS Workstation Graphics 
Software version 3.0, is used to draw and fill each panel 
on the monitor. The software may be run on other systems, 
with changes required in the screen plotting routines and 
the number of shades or color selection. Also the size of 
the input data files may be limited on other systems which 
do not have sufficient memory to store the input data in 
arrays. Results of this work show that the performance is 
limited by the plotting speed of the Micro-Vax. The 
calculation of the perspective image is approximately 11% 
of the total time required to complete the process. The 
remainder of time is used to draw the image on the screen. 

Using this software, an observer may select a position 
and heading and "see" the ocean bottom terrain in a 3-D 
perspective form. The observer heading is not limited, the 
position may be within the boundaries of the area 
selected, and the view is generated directly on the 
screen. Several options are available including a 
magnification factor, elevation scale exaggeration, saving 
the image to disk and the ability to respecify a new 
observer location. The speed of image generation is 
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directly affected by the resolution or size of the input 
data files, using 61 x 61 elevation points, an image is 
generated in 40 seconds. 

Chapter Two will cover the basic image geometry in 
order to develop the perspective projection transformation 
equation, and the orientation of the image plane. Chapter 
Three details the data file formats, resolution and how 
they are prepared for use in the main routine. Chapter 
Four covers the implementation of the algorithms in the 
program and any alternatives which were investigated. 
Chapter Five details program operation and Chapter Six 
presents conclusions and several recommendations for 
improvement or further study. 
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II. IMAGE FORMATION USING TRANSFORMATION GEOMETRY 



Inherent in this project is the ability to take the 
two dimensional gridded elevation data which represents a 
three dimensional object and display a perspective view on 
a two dimensional monitor. Two basic approaches, namely, 
parallel or perspective projection can accomplish this. 

A parallel projection is one where the points of the 
object are projected onto the image plane by parallel 
lines, that is the projection lines do not converge. This 
method has the advantage of maintaining relative 
dimensions of the object and is useful for obtaining 
measurements [Ref. 3, p. 52]. 

In the perspective projection all points are projected 
onto the image plane through a reference point which will 
be called the focal point [Ref. 4, p. 133]. In this 
method, the relative dimensions of the object are not 
preserved, but the image appears more realistic. Figure 
2-1 illustrates the two methods. Because the goal is 
feature recognition and interpretation through screen 
display and not precise measurement, the perspective 
projection is used. 
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Figure 2-1 Parallel and Perspective Projections 

[Ref. 3, p. 54] 
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The perspective projection may be generalized as a 
reference 3-D to 2-D transformation and is illustrated in 
matrix notation: 

A=[m]b (2.1) 

where: A-represents a vector in the image 3D 

coordinate system 

B-represents a vector In the object 3D 
coordinate system 

M-represents the transformation matrix 



A. COORDINATE SYSTEMS 

Equation 2.1 refers to the image 3-D coordinate system 
and the object 3-D coordinate system. The object is 
located in a right-handed 3-D cartesian coordinate system 
where each object point is described in terms of (X,Y,Z) 
geo-rectangular coordinates. The center of the earth is 
located at (0,0,0), the X-axis points to the intersection 
of 0 degrees latitude and 0 degrees longitude, the Y- 
axis points to the intersection of 0 degrees latitude and 
90 degrees east longitude and the Z-axis points to true 
north. Figure 2-2 illustrates this coordinate system. 
There is a direct relationship between the latitude, 
longitude, height with respect to sea level and the 
(X,Y,Z) object coordinates. The input data is supplied on 
a latitude, longitude grid and is converted to (X,Y,Z) 
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coordinates for use in the transformation equations. The 
image coordinate system is also a right-handed 3-D 
cartesian system. The focal point is located at 
(xO^yO^F). Figure 2-3 illustrates the image coordinate 
system. 



B. 3-D PERSPECTIVE TRANSFORMATION 

The elements of the [M] matrix in equation 2.1 are the 
cosines of the spatial angles that each of the image 
system axis x,y,z make with each of the object system axis 
X,Y,Z. This substitution results in the following specific 
form for the matrix [M] . [Ref. A, p. 139] 
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cosxX cosxY cosxZ 




^11 ”l2 ^13 


M = 


cosyX cosyY cosyZ 


= 


^21 ^22 ^23 




coszX coszY coszZ 




^31 ‘^32 ^33 



By convention the [M] matrix transforms from the 
object system to the image system [Ref. 4 , p. 139]. Using 
this matrix, every object point may be converted from the 
3-D XYZ coordinates to the xyz image coordinates. 

At this point the relation for a single object point 
is developed. Once defined, the relation may then be 
applied to all the vectors from the object to the focal 
point. Figure 2-4 illustrates the object and image space 
coordinate systems. Define the vector <FA> as a vector 
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OBJECT SPACE 




Figure 2-4 Object and Image Space Systems 



from the focal point (F) to a point on the image plane (A) 
and the vector <FB> as a vector from the focal point (F) 
to an object point (B) . Note that the image vector <FA> 
is collinear with the object vector <FB>. [Ref. 4 , p. 
141] 

FA = k FB 

where: FA - image vector 

FB - object vector (2.3) 

k - constant of proportionality 
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Using the image space coordinates (x,y,z) to describe <FA> 
and object space coordinates (X,Y,Z) to describe <FB> 
gives: 





*A - ^0 


FB *= 


■ Xb 


■ 


FA - 


va - V0 




■ ^F 




-f 




. ^B 


“ _ 



Then using equation 2.1 to express <FA> in the xyz image 
system yields: 

FA - k M j FB (2.5) 

or after substitution: 



^A - *0 




■ Xg - Xp ^ 


Va - V0 


■ k [ M ] 


Yb - Yf 


-f 




Zg - 2p 



Equation 2.6 is a form of the collinearity equation, it 
forms the basic relationship that the focal point, the 
image point and the object point are in a straight line. 
The values of xO and yO allow for the observer to be 
slightly misaligned with the image plane z-axis and 
generally are set equal to zero. This is the fundamental 
relationship used in this project to create the 
perspective views on the monitor screen. [Ref. 4, p. 141] 
This transformation will be called upon after the 
elevation points have been grouped into triangles and the 
observer location has been entered. Each elevation point 
in front of the observer will be transformed to image 
plane coordinates for subsequent drawing on the 2-D 
screen . 
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III. DATA FILE FORMAT FOR DIFFERENT SENSOR IMAGES 



Input to the routines consist of two data files, the 
bathymetric data and the image data. Each of these is 
described below. 

A. BATHYMETRIC DATA 

The bathymetric data was obtained from the Defense 
Mapping Agency/National Geophysical Data Center. The 
complete data base should be referenced as "Digital 
Bathymetric Data Base-Unclassified" or (DBDBU) . The area 
which was used here extends from 50 degrees north 
latitude, 140 west longitude to 32 degrees north, 120 west 
on a 5 minute by 5 minute grid. Each value describes the 
depth in meters below sea level at a given coordinate 
location. Values which lie above sea level are assigned 
-10 to avoid ambiguity. 

Extracting the area of interest is performed by a 
utility program, CROPELEV. FOR. The user enters the 
latitude and longitude of the southwest corner of the area 
of interest and the number of rows and columns desired. 
Each row and column is 5 minutes of latitude and longitude 
respectively. 
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After the area of interest is selected, the data is 
extracted from the original file and placed in row- 
column format where the rows correspond to latitude and 
the columns represent the longitude. Row numbers increase 
from south to north and column numbers increase from west 
to east. The first record of the file is a header 
containing the latitude and longitude (deg, min, sec) of 
the southwest corner, the latitude and longitude 
resolution in seconds and the number of rows and columns 
of data. Figure 3-1 illustrates this file format. This 
file is stored on disk and loaded into the variable 
RELEV(*,*) at runtime. 

B. IMAGERY DATA 

The sonar imagery data was obtained from the U.S. 
Geological Survey, Department of the Interior "Atlas of 
the Exclusive Economic Zone, Western Conterminous United 
States", [Ref. 5]. This atlas consists of 36 two degree 
mosaic images at a scale of 1:500,000. Coverage extends 
from 49 degrees north, 130 west to 30 degrees north, 117 
west. These images were obtained using a unique side-scan 
sonar system called GLORIA (Geological Long-Range Inclined 
Asdic). [Ref. 5, p. 2] This sonar system allows mapping 
of large areas of ocean bottom on a single pass of the 
ship. Figure 3-2 lists several characteristics of the 
GLORIA system. 
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Figure 3-1 Bathymetric Data File Format 
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Gloria Side Scan Sonar Specifications 



Size 7.75 X .66 meters 

Weight 2 tons 

Scan 30 km/side = 60 km 

swathwidth at 10 kts. 

Power 10 kw/side 

Beamwidth 2.5 deg azimuth 

10 deg vertical 

Resolution 30 x 218 meters/pixel 

at 5000 meters depth 

Frequency 6.5 kHz 



Figure 3-2 Some Characteristics of the 
Gloria II Side-Scan Sonar 
[Ref. 6 , p. 7] 



Each mosaic image is a half-tone black and white 
print of the acoustic reflectance of the sea floor with 
white representing the highest reflectance and black the 
lowest reflectance. As seen in Figure 3-2, the resolution 
of each pixel is distorted (30 meters by 218 meters) due 
to the ships motion along a track perpendicular to the 
scanning direction. The smaller value (30 meters) is the 
resolution in the cross-track direction while the larger 
value (218 meters) is in the along track direction. Prior 
to forming the mosaic images, aspect ratio distortion is 
removed by correcting for geometric and radiometric 

17 
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distortions in the raw data. The final resolution, after 
correction, used to create the mosaics is about 50 meters 
X 50 meters. This relatively low resolution image can 
show only large scale features but may point out areas 
which deserve additional attention. Also note that this 
resolution is several orders of magnitude better than the 
available bathymetric data. [Ref. 6, p. 3] 

Digital sonar imagery data was not available for use 
in this thesis. In order to determine the effectiveness of 
this imaging method, the digital image data was created by 
locally digitizing the mosaics contained in the atlas. 
This provided digital imagery data in the following form, 
512 X 512 pixels with each pixel using one byte of storage 
resulting in 256 shades of gray. This data is stored on 
disk in direct access form as 512 fixed length records 
with 512 bytes per record. Record one is the northern end 
of the image and record 512 is the most southern end of 
the image. Figure 3-3 illustrates the orientation of the 
image file. The disk file is loaded into the variable 
IMAGE ( * , * ) at runt ime . 

Ideally, the digital imagery data would be directly 
available, either from the GLORIA side-scan sonar or from 
other sonar systems. If this were the case, the data would 
require processing to remove the different sources of 
error and to convert it from its native form into the 
form required in this project. Processing techniques for 
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Figure 3-3 Image File Orientation 



GLORIA digital data is discussed in Ref. 1 , "Processing 
Techniques for Digital Sonar Images from GLORIA" by Pat S. 
Chavez, Jr. 

The current limitation for input data is 70 rows by 70 
columns for the elevation data and 1024 rows by 1024 
columns for the image data. These dimensions may be 
increased to accommodate higher resolution data; the 
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tradeoff is speed of execution. The number of triangular 
panels is a function , of the elevation data size, 
specifically: 

#panels = ((#rows - 1)* 2) * (#cols - 1) (3.1) 

An elevation file with resolution 20 x 20 has 722 
panels; increasing the resolution to 100 x 100 results in 
19,602 panels. Each panel is projected point by point then 
sorted and drawn for every image view constructed. The 
increase in computation time follows directly. In addition 
to the elevation file resolution, the image file 
resolution may be increased. The cost in computation is 
not severe in this case because the image file is accessed 
only once when the image pixels are averaged to calculate 
the color or gray level for the individual panels. This is 
done prior to displaying the first image and is not 
repeated unless the program is exited and restarted. 
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IV. ALGORITHMS FOR COMBINED 3-D PERSPECTIVE DISPLAY 



The main routine begins by reading two input data 
files into the arrays RELEV(*,*) and IMAGE(*,*). Figure 
4-1 illustrates the steps involved in displaying the image 
on the monitor. 

The procedure begins by forming the triangles in the 
grid of the elevation data, averaging the image pixels 
inside each triangle and assigning the average gray 
intensity to the panel (one triangle) . The observer 
location is read, then the image plane coordinate system 
is constructed. Based on the orientation of the image 
plane the [M] matrix is calculated and used to map the 
elevation points to the image plane. A 2-D to 2-D affine 
transform is used to convert the points from image plane 
coordinates to the screen coordinates. To display the 
perspective view on a flat screen, triangles which are 
hidden behind foreground objects must be removed. This 
hidden surface removal is performed by calculating the 
distance from the panel to the observer and drawing the 
panels in sequence from the farthest to the nearest. This 
chapter will cover these areas in detail and also discuss 
alternatives which were considered, but not implemented. 
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Begin 



Load input data file 

Read in elevation (bathyroetric) data 
Read in imagery data 

Construct triangles in elevation data grid 
Average the gray shade in each triangle 

Get observer location 

Form the image plane vectors 

Calculate the elements of the [M] matrix 

Transform to the image plane 

Determine which panels are visible 

Convert to screen coordinates 

Remove hidden surfaces 

Plot results on screen 

End 



Figure 4-1 Basic Program Flow 



A. POLYGON FORMATION AND SHADING 

Solid objects may be represented in numerous ways. 
Some objects lend themselves to being described in terms 
of a number of planes or surfaces, A cube, for example may 
be precisely defined with six planes [Ref, 8, p, 189], As 
the object or scene becomes more complex, the number of 
surfaces required to accurately describe it increases. 
This method of representing a 3-D surface by plane 
surfaces lends itself particularly well to this 
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application. It allows for easy calculation of the color 
or gray level by averaging of pixels contained within the 
panel, easy transformation to the image plane, and easy 
use of graphics hardware to execute the polygon fill 
operation. This last feature is most important; earlier 
work in this area pointed out the amount of time required 
to draw based on a point by point method [Ref. 2, p. 56]. 
Using the graphics hardware polygon fill capability 
alleviates this problem. 

The input elevation data is supplied in gridded form 
and is easily partitioned into squares and ultimately 
into triangles. Figure 4-2 illustrates the partitioning of 
the elevation data into triangles. The procedure for 
selecting the gridsquares and forming the triangles is 
given as psuedocode in Figure 4-3. The column coordinates 
of each point are stored in IA(*) and the row coordinates 
are stored in the JA(*) array. Starting in the southwest 
corner, the program selects four elevation points which 
form a gridsquare. These four points are called node_a, 
node_b, node_c and node_d. The gridsquare is divided into 
two triangles by calculating the equation of the line 
which divides it. This dividing line is completely 
described by the slope and y-intercept. The column 
boundaries are defined by IA(node_a) and IA(node_b) ; the 
row boundaries are defined by JA(node_a) and JA(node_c) . 
For each incremental step from the eastern (right hand) 
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View from above looking down on the terrain. 

-Terrain elevation points are connected 
to form triangular polygons with common 
edges . 



Figure 4-2 Polygonal Terrain Construction 

[Ref. 3, p. 35] 
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BEGIN 



DO select a gridsquare 

(*select four nodes which make one gridsquare*) 

READ node_a 

READ node_b 

READ node_c 

READ node_d 

(*calc slope and y-intercept of line forming 
the triangle*) 

slope = rise / run 

y-intercept = y - slope * (x) 

DO M = IA(node_b) to IA(node_a) step -1 

IY= (slope) * M + y-intercept 

DO L = JA(node_b) to lY step -1 

(*average image pixels in lower triangle*) 

END DO 

DO L = lY to JA(node_c) step -1 

(*average image pixels in upper triangle*) 
END DO 
END DO 

PL_AR(*,*) = nodes of triangle, average gray shade 
END DO 
END 

Figure 4-3 Polygon Formation Psuedocode 
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boundary, IA(node_b) to the western (left hand) boundary, 
IA(node_a) , a value of y equal to lY is calculated based 
on the equation of the dividing line. The image pixels 
above and below this line are averaged separately for the 
upper and lower triangles. The process is repeated for 
each incremental step in the columns from right to left. 
This scanning and averaging pattern is illustrated in 
Figure 4-4. Note that the scan pattern is from bottom to 
top, right to left in each triangle. [Ref. 2, p. 39] The 
average shade of the triangle along with the panel number 
and vertices are stored in array PL_AR(*,*). The data 
structure is shown in Figure 4-5. Notice that PL_AR(*,5) 
generally is not used, it will be used to store the 
maximum distance of the plane from the observer. 

B. IMAGE PLANE FORMATION, PERSPECTIVE VIEW CALCULATION 

After the shading for each panel is calculated, the 
observer's location is entered. This consists of latitude 
and longitude (deg, min, sec) , height with respect to sea 
level and heading. The latitude, longitude and height data 
is converted to object space coordinates X1,Y1,Z1 while 
the heading is used to select a second point in front of 
the observer. This second point X2,Y2,Z2 forms the line of 
sight (LOS) vector for the observer. The negative LOS 
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Figure 4-4 Gridsquare 
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TRIANGLE IDENTIFIER 


NOOE_A 

FIRST VERTEX I.O. 


NOOE_0 

SECOND VERTEX I.O. , 


NOOEj: 

THIRD VERTEX I.O. 


CRAY SHADE 


MAXIMUM DEPTH 


1 


3660 


3721 


3720 


142 


109979 


2 


3659 


3720 


3719 


122 


109178 


3 


3659 


3660 


3720 


123 


109178 


4 

• 


3599 

• 


3660 

• 


3659 

• 


117 

• 


108434 

• 


• 

• 


• 

• 


• 

• 


• 

• 


• 

• 


• 

• 



Figure 4-5 Layout of PL_AR Array 



vector is the z-axis of the image coordinate system as 
seen in Figure 2-3. This is formulated as shown; 



ZVECj^ = XI - X2 
ZVECy = Y1 - Y2 
ZVECg =21-22 



2VEC = 2VECx X + 2VECy 



(4.1) 

y + 2 VEC 2 z 
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The image plane y-axis is constructed by creating a vector 
which points upward from the center of the earth to the 
observer location. The (X,Y,Z) coordinates of the observer 
location (X1,Y1,Z1) form this vector. The y-axis is 
described by: 



YVECx 


= 


XI 




YVECy 


= 


Y1 


(4.2) 


YVECg 


= 


Z1 




YVEC = YVECx 


A. 

X 


+ YVECy y + YVEC 2 z 





The image coordinate system is a right-hand system. 

Therefore, the cross product of the y-axis and z-axis 

forms the x-axis. The x-axis of the image plane is 

constructed as follows: 

XVECx = ((YVECy x ZVECz) “ (YVECz x ZVECy)) 

XVECy = ((YVECg X ZVEC^) - ( Y VEC^ x ZVECg)) (4.3) 

XVECg = ((YVECx X ZVECy) - (YVECy x ZVEC^ ) ) 

XVEC = XVECx X + XVECy y + XVECg z 

The [M] matrix is calculated using the relationship 
developed in Chapter Two, equation 2.2. Recall that each 
element of the [M] matrix is the cosine of the spatial 
angle between the axis of the image coordinate system and 
the object coordinate system axis. The object space 
coordinate axis may be represented as unit vectors: 

XAXIS = lx + 0y + 0z 

YAXIS -= 0x + ly + 0z (4.4) 

ZAXIS = 0x + 0y + Iz 
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The cosine of the spatial angle between two vectors is 
determined from the dot product of the two vectors. 



A • B 



B cos & 



AB 



The dot product may be expanded as shown; 



cos © 



AB 






lA i IB 



cos 6 



AB 



lA I IB I 



(4.5) 



(4.6) 



Substituting the image space x-axis vector (xvec) and the 
object space X-axis vector (XAXIS) for [A] and [B] allows 
Mil to be calculated: 

cos e.v « XVECx-XAXISx ^ XVECyXAXISy XVEC;, -XAXIg^ 

IXVECI IXAXISI 



(XVECy-1) + (XVECy‘0) + (XVECg*®) 
IXVECI • 1 



(4.7) 



H 
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XVECy 
IXVEC I 



The other elements of the [M] matrix may be found in a 



similar fashion. This results in: 



Mil - 



M31 



XVECy 
IXVEC I 

YVECy 
lYVEC I 

ZVECy 
IZVEC I 



Mi2 = 

«22 - 

H32 " 



XVECy 
IXVEC I 

YVECy 
lYVEC I 

ZVECy 
IZVEC I 



Mi3 = 

«23 " 
H33 ■ 



XVEC2 
IXVEC I 

YVECg 
lYVEC I 

ZVECz 
IZVEC I 



(4.8) 
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This method of calculating the image coordinate axis 
and the corresponding [M] matrix elements results in a 
viewing plane which is oriented vertically and faces in 
the direction of the observer course. 

An alternative to this method is that, instead of 
specifying the course to see the perspective view, one can 
have the image plane always face the object. This would 
simulate rotating the object while keeping the observer 
location fixed. The decision to use the first method is 
based on the feeling that it results in a more natural 
view for the observer. 

C. TRANSFORM TO THE IMAGE PLANE 

After the [M] matrix is constructed the elevation 
points are projected onto the image plane. Points which 
lie behind the image plane or behind the observer are not 
projected. The method used to select which points to 
project and the projection equations are presented next. 

Determining which elevation points are located behind 
the image plane is done by finding the cosine of the 
angle between a vector formed by the negative z-axis of 
the image plane and a vector drawn from X2,Y2,Z2 to the 
object point. Solving for the cosine of the angle and not 
the angle itself eliminates the need to use trigonometric 
functions, and provides faster execution. Psuedocode for 
this procedure is given in Figure 4-6. 
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BEGIN 



DO IR = 1 to Number of elevation pts. 

(*get XYZ coordinates of point*) 

X = XYZ(IR,1) 

Y = XYZ(IR,2) 

Z = XYZ(IR,3) 

(*form the observer LOS vector*) 

OBS_VECX = X2 - XI 
OBS_VECY = Y2 - Y1 
OBS_VECZ = Z2 - Z1 

(*form the object vector*) 

OBJ_VECX = X - X2 
OBJ_VECY = Y - Y2 
OBJ_VECZ = Z - Z2 

(*find the cosine of the angle between the 
OBS_VEC and OBJ_VEC*) 

COS_THETA = (dot product of OBS_VEC, OBJ_VEC)/ 

(Mag(OBS_VEC) x Mag (OBJ_VEC) ) 

IF (COS_THETA > 0) THEN 

Project onto image plane 
Calculate depth of panel 

ELSE 

Mark as a non-visible point 
END IF 

END DO 



END 



Figure 4-6 Psuedocode Determine Which Panels to Project 
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The procedure begins by forming two vectors, the 
observer LOS vector and the observer to object vector 
(object vector) . Ideally, the origin of the image plane 
coordinate system would be used to form the "observer to 
object vector", but the (X,Y,Z) coordinates of this point 
are not readily known. The observer location (X1,Y1,Z1) is 
not used because it is behind the image plane and will 
result in some panels being viewed which are behind the 
image plane. The point (X2,Y2,Z2) is used as an 
approximation of the origin of the image plane. If the 
angle between the LOS vector and the object vector is 
greater than 90 degrees, the point lies behind the image 
plane and will not be projected onto the image plane. 
Figure 4-7 illustrates the orientation of the vectors and 
the image plane. The cosine of the angle between the two 
vectors is again found by using the dot product method. 

If the angle is less than 90 degrees, the cosine will 
be a positive number. If the angle is greater than 90 
degrees the cosine will be negative; this indicates the 
point is not in the hemisphere in front of the observer. 
If the point is not visible in the forward hemisphere then 
the point is flagged as hidden. Generally each point is 
associated with up to six triangles; marking each point 
as hidden will prevent attempting to draw triangles which 
may be partially behind the image plane. 
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Figure 4-7 Image Plane Object Vector 



The points which are in the forward hemisphere are 
projected onto the image plane using the relation 
developed in Chapter Two, equation 2.6 which is repeated 
here for clarity; 



- ^0 




■ Xb - Xp ■ 


va - V0 


= k [ M ] 




-f 




Zb - Zp 
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Expanding to three separate equations gives: 



*A - *0 


- k[ 


^1 1 ( Xg Xp ) 




^13^2g-2p) J 


Va - V0 


- k[ 


M2i(Xg-Xp) 




^23^^B"2f^ ] 


-f 




Mai^Xfi-Xp) 


+ M32(Yb-Yf> 


^33^^B”2f^ ] 


Then dividing (4.10) and 


(4.11) by (4 


.12) yields: 


V _ V 


m — f 


MilCKg-Xp) 




+ Mj^3(2g-2p) 


0 




M32^(Xg-Xp) 




+ M33(2g-2p) 




B — f 


M2i(XB-Xp) 


+ M22(Yb-Yf> 


M23(2B“2p) 


y ^0 


B — ^ 


M3i(XB“Xp) 


M32(Yg-Yp) 


M33(2B“2p) 



( 4 . 13 ) 



( 4 . 14 ) 



Equations (4.13) and (4.14) are the relations used to 
transform the elevation data points to the image plane. 
They allow projecting all elevation points in the 
hemisphere forward of the observer at one time. Altering 
the size of the image plane will not require the elevation 
points to be projected a second time. [Ref. 4, p. 142] 



D. AFFINE IMAGE PLANE TO SCREEN TRANSFORM 

The image plane coordinates must be transformed to 
the screen coordinate system. Figure 4-8 illustrates the 
image plane and the screen coordinate systems. The maximum 
frame size values of the x and y axis in the image plane 
system are determined by the desired magnification or 
field of view and the y-scale factors. 

In general, the affine transform will map between any 
2-D coordinate systems allowing for a rotation of the 
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Figure 4-8 Image to Screen Coordinates 



axis, horizontal and vertical scale changes, two 
translations and a non-perpendicularity of the axis. The 
general form of the affine transform is [Ref. 4, p. 593]: 



*1 


1 


^2 


-bl ■ 




^2“ ^1 


. yi . 


^1^2 “ ^2^1 


. ~^2 


®1 . 




. y2~ ^2 . 



where: (cos ff -i sin P) 

= Sy (-sin P) 

C represents a non-perpendicularity of the axis 
^ represents a rotation of the axis 
= X translation of origin 
C2 = y translation of origin 
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starting with a magnification factor of 1.0 and a y- 
scale factor of 1 . 2 , the x and y frame sizes are 35000 
and 29000 respectively. These values must be mapped onto a 
screen coordinate system which ranges from 1-512 in both 
the X and y direction. There are two scale changes, two 
translations, no rotations and no non-perpendicularities 
involved with this transformation. This results in sigma 
equal to zero and beta equal to zero. Substituting these 
values results in: 



■ s. 



*2 ' ® “2 - “y 

Inserting these values and expanding the matrix equations 



0 

S, 



( 4 . 16 ) 



gives: 



- C- 






ll. 



- c. 



( 4 . 17 ) 



The scale factors are formed at runtime to allow the 
changes in magnification (field of view) and elevation 
scale. Figure 4-9 is psuedocode which implements the 2-D 
transform. Note that the variables Al, A2 , Bl, B2 are in 
terms of the image plane frame size. This approach allows 
the user to change the magnification or elevation scaling 
interactively at this point. Also note that each elevation 
point is checked for non-visible status prior to the 
transform. This checking saves execution time by avoiding 
points not in the field of view. 
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BEGIN 



Cl = image x axis maximum 

C2 = image y axis maximum 

A1 = image x axis frame size / screen x axis maximum 

B2 = image y axis frame size / screen y axis maximum 

DO IR = 1 to Number of elevation pts, 

IF (point .EQ. visible) THEN 
xscreen = (ximag - Cl) / A1 

yscreen = (yimag - C2) / B2 

END IF 

END DO 

END 



Figure 4-9 Affine Transform Psuedocode 



E. HIDDEN SURFACE REMOVAL 

In order to properly display the perspective view, the 
surfaces hidden behind front surfaces must be dealt with. 
Earlier work utilized a Z-buffer to perform the hidden 
surface removal [Ref. 2 , p. 41]. The Z-buffer method 
utilizes a pixel by pixel depth buffer. Each pixel 
position on the screen corresponds to a location in the 
buffer. The buffer is initialized to a maximum depth, then 
prior to drawing each pixel the depth is compared with the 
depth already stored in the buffer. If the pixel depth is 
less than the depth stored in the buffer, then the pixel 
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is drawn and the new depth is stored in the depth buffer. 
Earlier work by Coleman pointed out the poor performance 
of this technique when executed in software [Ref. 2, p. 
55] . Several other methods of hidden surface removal were 
investigated; each will be discussed briefly and the 
disadvantages listed. 

First, the method of "bounding rectangles" was 
investigated. In this method a rectangle is formed around 
each triangle. This rectangle is as small as possible and 
is constructed from the three vertices forming the 
triangle. The (x,y) coordinates of each vertex are checked 
and the minimum and maximums of both the x and y 
coordinates form the boundaries of the rectangle. Now the 
vertices of all other triangles are compared to the 
rectangle. If the (x,y) coordinates of a vertex lie within 
the boundary of the rectangle, the triangle is marked as 
overlapping and the Z-buffer routine is used to handle the 
hidden surface removal. If the triangle does not overlap 
any other triangles, it is drawn to the screen. Originally 
the intent was to reduce the number of panels which needed 
to be drawn using the Z-buffer method. This was not 
achieved because the panels are triangular vice 
rectangular. Several experiments showed that no triangles 
passed the bounding rectangle test and consequently no 
savings in execution time was realized. [Ref. 9, p. 246] 
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A second method utilizes a test function for each side 
of a triangle. A test function is formed for each side of 
a triangle. This test function is formed from the equation 
of the line describing the side. The vertices of all other 
triangles are checked against the test function. This 
results in nine comparisons for every triangle checked. 
The calculation time is excessive, also the results were 
similar to the results of the bounding box test. Very few 
triangles passed the test and the Z-buffer method had to 
be used on the majority of panels, with no savings in 
execution time. [Ref. 9, p. 24] 

Another approach, known as the "Painter's algorithm", 
is not hidden surface removal at all. [Ref. 8, p. 265] 
Here the panels are drawn to the screen in sequence from 
background to foreground. The panels which are hidden from 
view are covered with the panels which are closer to the 
observer. This is the approach used in this software. Each 
panel is sorted into order based on the maximum distance 
of the vertices from the observer location. The choice of 
sorting algorithms has a large affect on the efficiency of 
this method. Originally, the routine incorporated a simple 
bubble sort. The cost in performance is not noticeable 
with a small number of panels, but as the resolution of 
the bathymetric data increases, the number of panels 
increases as shown in equation 3.1 which is repeated here. 
# panels = ((#rows - 1)*2) * (#cols - 1) (4.18) 
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The poor performance of the bubble sort becomes apparent 
as the number of panels increases beyond 500. Two 
alternate sorting routines were used with this routine, 
the shell sort and the quick sort. Table 4-1 is a table 
which shows the measured performance of these sort 
routines. Clearly for large numbers of panels the quick 
sort is the most efficient. Based on these results the 
quick sort was selected for use in the hidden surface 
routine. Figure 4-11 is psuedocode which describes the 
hidden surface removal procedure. 



TABLE 


4-1 SORT 


PERFORMANCE 






Time 


in seconds 




# panels 


Bubble 


Shell 


Quick 


500 


2 


<1 


<1 


1000 


8 


<1 


<1 


2000 


33 


<1 


<1 


3000 


80 


1 


1 


5000 


210 


2 


1 


10000 


900 


5 


1.5 


20000 


— 


14 


3 


50000 


— 


56 


9 
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BEGIN 



(* create sort key by counting the panels in front of 
the image plane. Depth_key(L) contains the panel ID 
number * ) 

DO IR = 1 to NPLANES 

IF ( panel .NOT. behind field of view) THEN 
INCREMENT COUNT 
Depth_key( COUNT) = IR 
END IF 

END DO 

(* calc maximum depth of each visible plane *) 

DO L = 1 to COUNT 
IR = Depth_key(L) 

PL_AR(IR,5)= 

MAX (DEPTH (node_a) , DEPTH (node_b) , DEPTH (node_c)) 



END DO 

(* call sort routine to sort the panels on the 
maximum depth*) 

CALL QUICKSORT 



END 



Figure 4-11 Psuedocode for Hidden Surface Removal 
by Sorting 



Upon completion of the hidden surface removal routine, 
the screen is erased and the panels are plotted. The 
vertices for each panel along with the shading are stored 
in array PL_AR(*,*). The shading value is clipped to a 
range of 1-250 vice the original 1-256 due to the DEC VMS 
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window manager. This maintains the background color of the 
screen and the window attributes. The panels are read in 
sequence using the sort index. The (x,y) screen 
coordinates of each vertex and the shade are passed to the 
plot subroutine and the panel is drawn on the screen. 



F. RESULTS 

Data cropped out of the primary bathymetry data file 
covers an area bounded by 37:00:00 N. latitude, 126:00:00 
W. longitude to 36:00:00 N. latitude, 125:00:00 W. 
longitude. Figure 4-12 is the digitized image of the 
mosaic used for the image data. The first synthesized view 
Figure 4-13, is from an observer location of 36:0:0 N. 
125:30:0 W. , 4000 meters below sea level and heading 
equal to 000 degrees. The y-scale exaggeration is six and 
the magnification factor is one. Figures 4-14 a to c are 
displays of the same area from different observer 
locations. Figures 4-14a and 4-14b are offset from the 
original location east and west by 25 minutes of 
longitude. Figure 4-14c is from the same location as the 
original but the height is 2000 meters below sea level. 
These images display the ability to "see" the ocean bottom 
in a new fashion. The resolution of the input elevation 
data was originally 5 minutes by 5 minutes. This data was 
interpolated in both latitude and longitude to a 
resolution of 1 minute by 1 minute in order to achieve a 
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better shading rendition. This reduces the number of 
pixels averaged for each triangle and improves the shading 
of the generated image. Figure 4-15 is an image using the 
original 5 minute by 5 minute resolution data. Notice 
that the size of the panels is larger and the shading 
rendition is much more coarse. In addition to affecting 
the shading, the resolution of the elevation data directly 
affects the physical shaping of objects in the terrain. 
This program does not utilize any method of curve fitting 
between adjacent elevation points and simply draws 
straight lines between the points. This results in 
objects being coarsely approximated in the synthesized 
view. 
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Figure 4-12 Mosaic Image 



Figure 4-13 



Obs Loc: 36:00:00 N Depth 4000 m 
125:30:00 W Heading 000 deg 
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Figure 4-14a Obs Loc 36:00:00 N 

125:05:00 W 



Depth 4000 m 
Heading 340 deg 




Figure 4-14b Obs Loc 36:00:00 N Depth 4000 m 

125:55:00 W Heading 020 deg 
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Figure 4-15 5x5 Minute Resolution 



Figure 4-14c Obs 



Log 



36:00:00 N Depth 2000 m 
125:30:00 W Heading 000 deg 






V. OPERATIONAL ASPECTS OF THE PROGRAM 



The software utilizes an input file to load the names 
and geographic locations of the data files. This input 
file may be created by an ASCII editor or created within 
the program by following the prompts. The format of the 
input file is shown in Figure 5-1. 

If the input file exists, the user simply enters the 
name when prompted. The program will seem to pause at this 
point while reading in the data files, forming the panels 
and averaging the image pixels. When complete with these 
steps, the software will prompt for the observer location. 
This location is entered in terms of observer latitude and 
longitude in degrees, minutes and seconds. Southern 
latitudes and western longitudes are entered as negative 
values. Also entered is the elevation of the observer and 
the heading of the observer. Elevation is in meters with 
respect to sea level with values greater than sea level 
entered as positive values and elevations below sea level 
entered as negative values. The heading is entered as a 
value from 0-360 degrees where 0 deg. is north, 90 deg. is 
east, 180 deg. is south, 270 deg. is west and 360 deg. is 
north again. 
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Elevation data file name 

# of rows of elevation data 

# of columns of elevation data 
Image data file name 

# of rows of image data 

# of columns of image data 

Latitude/Longitude of Northwest corner of image 
data in (deg, min, sec) 

Latitude/Longitude of Southeast corner of image 
data in (deg, min, sec) 

Title for the screen image 

Figure 5-1 Input File Format 



At this point the software will draw the perspective 
view on the screen. It takes 40 seconds to produce one 
image using 61 x 61 elevation points and 512 x 512 image 
points. After the image is drawn a menu is presented. The 
operator may alter the magnification or field of view, 
change the elevation scale exaggeration, enter a new 
observer location, save the image or exit the program. 

A. ALTER MAGNIFICATION OR FIELD OF VIEW 

The perspective view transformation simulates viewing 
through a camera viewfinder. The observer selects the 
location and points the camera, and the image is formed in 
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the viewfinder. The magnification or field of view is 
changed by changing the focal length of the lens. 

The magnification factor of the initial image is one. 
The values of the image plane x-axis and y-axis maximums 
are initially set to 35000 and 29000 respectively. This 
loosely corresponds to a 35mm x 29mm frame size of a 35mm 
camera. The focal length corresponds to an initial focal 
length of 50mm. The field of view is given by: 

f xframe 

Field of View « 2 tan“^ L2( focus) J 

Using the values above results in a field of view of 36.8 
degrees. This field of view is approximately equal to the 
field of view from a 35mm camera with a 50mm "normal" 
lens. This field of view corresponds to a magnification 
factor of one. Table 5-1 shows the field of view for 
various magnifications. 

When changing the magnification on a camera, the 
standard method is to change the focal length of the lens. 
This method presents a disadvantage here in that it 
requires the complete 3-D to 3-D transform to be 
performed for every magnification change. By changing the 
frame size instead of the focal length, the new image is 
generated faster because the change is made during the 
affine image plane to screen plane transformation. The 
magnification changes are entered as values with respect 
to the normal magnification. 
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TABLE 5-1 MAGNIFICATION VERSUS FIELD OF VIEW 



Magnification 



Field of View (deg) 



. 5 

1.0 

1.5 
2 . 0 

2 .5 

3 . 0 

3 . 5 

4 . 0 

4 . 5 

5.0 



69.9 

38.5 

26.2 

19.8 

15.9 

13.3 

11.4 

10.0 

8.9 

8.0 



B. CHANGING THE ELEVATION SCALING 

Altering the y-axis scale allows exaggeration or 
reduction of the y-scale. Views of the terrain with the 
elevation scaling set at 1:1 result in views which have 
very little appearance of height. By modeling the frame 
size as a 35mm camera frame, a 1.2:1 scaling is 
introduced. Research conducted by the U.S. Army Research 
Institute for the Behavioral and Social Sciences indicated 
that vertical scaling ranging from 1.25:1 to 1.50:1 
presented the most natural views. [Ref. 3, p. 37] 

The elevation scaling is changed in a fashion similar 
to the magnification as a value with respect to the 1:1 
scale. 
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C. ENTER A NEW OBSERVER LOCATION 



This option allows the operator to enter a new 
observer location. The current settings for magnification 
and elevation scaling are maintained. 

D. SAVING AN IMAGE 

Selecting this option allows multiple images to be 
saved to a disk file. When the option is selected the 
first time, it will prompt for the file name for the 
images to be stored in. Two data files are created, *.IMG 
and *.SIZ. If the filename extension is specified then it 
is used instead of the .IMG extension. Images will be 
saved to *.IMG file for the entire session each time the 
(S)ave option is selected. The entire set of images may be 
viewed sequentially using a program named PLAYBACK. FOR. 
This will display all the images stored in the image file 
at a rate of one image per second. Another program, 
UIST0512.F0R will convert a single frame from the saved 
images to a 512 x 512 row-column format, with 512 fixed 
length records and 512 bytes per record. The *.SIZ file 
contains the buffer size data required by the playback 
routine. 
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VI. CONCLUSIONS 



A. GENERAL 

The goal of this thesis was creating the 3-D 
perspective view of ocean bottom terrain and applying 
shading or color information from an image file. This goal 
was achieved and the ability to "see" the terrain in this 
fashion holds promise for future work in many fields. 

The program presented here combines terrain elevation 
data with shading information from a second data source. 
This fusion of bathymetry and imagery data results in a 
3-D perspective view of the ocean bottom in a fashion 
which appears natural to the observer. This type of 
display system, which combines multiple sources of data 
into a form that is easier to comprehend, has potential 
uses in mapping, geologic exploration and weapon system 
displays. It has direct use in any geographic information 
system where the data stored consists of (x,y) coordinates 
and an attribute. This program could merge any two layers 
of information from such a system and display a view which 
is easier to interpret. The data utilized in this project 
was not a part of such a system, but easily could have 
been. In this case the first layer of data would be the 
(x,y) location and the elevation with respect to sea level 
and the second layer of data would be the (x,y) location 
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and the shading value. Combining the layers of data from 
multiple sources is a form of sensor fusion. 



B. LIMITATIONS 

There are several limitations in this program. First, 
the image plane orientation is fixed in the vertical 
direction with respect to the observer heading and 
location. This implies that the perspective view formed 
simulates the view from a vehicle in level flight. 
Allowing the image plane to deviate from vertical would 
allow the observer to look down on the terrain and 
generate views not presently permitted. Second, the 
shading values are clipped to 250 vice 255 shades. This is 
done in order to maintain five separate entries in the 
virtual color map which correspond to the system colors. A 
third problem exists with the geographic areas of the 
input data. The conversion from latitude/longitude 
coordinates to geo-rectangular coordinates is not 
implemented at the hemisphere boundaries such as at the 
equator (0 deg. N.), or at 0 deg. longitude. The input 
data must be completely above or below the equator and 
similarly completely east or west of 0 deg. longitude. 



C . PERFORMANCE 

Forming the 3-D perspective views requires the input 
data to be transformed from one 3-D coordinate system to 
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the image plane system and ultimately to the screen 
coordinate system. As discussed in Chapter Three, the size 
of the input data files, specifically the elevation data 
will affect the performance of the program. Table 6-1 is 
a breakdown of the amount of time required to complete 
each of the major sections of the program. The overall 
performance is satisfactory and breaks down to 40 seconds 
to complete a view after the observer location is entered. 
The plotting of the view on the screen is 89% of this 
time, while the calculation involved accounts for only 11% 
of the time. Clearly the plotting time must be improved 
prior to spending an appreciable amount of time in the 
calculation subroutines. The length of time required by 
the plotting routine is related to the number of triangles 
being passed to the graphics hardware. The communication 
involved in passing the vertices of each panel to the 
graphics processor is the bottleneck. Improvements may be 
achieved in this area by writing a driver to directly load 
the appropriate values to the graphics hardware and 
bypassing the overhead of the software routines presently 
utilized. 

The combination of the perspective transformation and 
the hidden surface removal routines account for 95% of the 
11% calculation time as seen in Table 6-1. Other methods 
of achieving slight performance increases may include a 
hardware implementation of the perspective transformation 



55 



TABLE 6-1 COMPUTATION 


TIME 


Subsection 


Resolution 




61 X 61 


13 X 13 




(time 


in seconds) 


Set up graphics 


.80 


.53 


User input 


— 


— 


Read elev. file 


.25 


. 11 


Read image file 


10.22 


10.69 


Form triangles 


3.3 


. 3 


Average gray shade 


6.87 


2.54 


Observer location 


— 


— 


Form [M] matrix 


.01 


. 01 


[M] multiply 


2.51 


. 12 


Affine transform 


.23 


. 01 


Hidden surface removal 1.79 


. 07 


Plotting 


35.24 


1.76 


Total elapsed after 
Observer location 


39.78 


1.97 


Calculation time 
percentage 


11.4% 


10.6% 



matrix multiplies and a more efficient hidden surface 
removal routine. Altering the hidden surface removal 
routine most likely will result in more calculation time 
required in that area, but may result in substantial 
savings in the plotting time. The approach for hidden 
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surface removal used here does not actually remove hidden 
triangles but draws them in order from the background to 
foreground. By altering the hidden surface routine so that 
hidden panels are not plotted at all, a savings may be 
made in the time required to plot the view. 



D. FUTURE WORK 

This program is a first step in the direction of 
combining various layers of data into a single form. There 
are several areas which are being considered for further 
work. 

The type of data merged onto the elevation grid is not 
limited to the imagery data used here. Magnetic anomaly 
data is available and could also be merged onto the grid. 
This would permit the observer to make correlations 
between the terrain and the magnetic response. 

The original effort in this area began with aerial 
photographs and surface terrain. Applying the terrain 
elevation data and an aerial photograph to this program 
would broaden the versatility of the program greatly. 

Combining two layers of data is a starting point, the 
ability to combine three or four layers onto the 
elevation grid using colors to distinguish the layers may 
be possible. 

In order to make this program truly useful, a user 
interface must be designed. As a minimum this requires a 
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method of entering observer locations based on the current 
location and perspective view, A particularly effective 
user interface would involve using a pointer to select an 
area on the screen and tell the program "I wish to go 
there” . This would permit "driving around" an area of 
interest . 

Finally, a method of loading data from the storage 
device and merging it with the currently loaded data is 
required. This would remove the boundaries from the image 
displayed on the screen so that when the observer location 
is near the edge of the current area, the adjacent areas 
will be loaded and displayed. 
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APPENDIX A. PROGRAM SUMMARY 



This appendix lists the different routines of the 
program by subroutine name and briefly lists the inputs 
and outputs of the subroutine. Also given are the names of 
the subroutines which interact with each subroutine. 

1. S0NAR3D 

A. Function 

Main routine, calls other supporting routines to 
complete the 3-D perspective view 

B. Input 
None 

C . Output 
None 

D. Calling routine 
None 

E. Called routine 
INPUT 
READ_ELEV 
READ_IMAGE 
TER_XYZ 
IM_REFAVG 
OBS_LOC 

M ORIEN 
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NEW IJ 



XY2IJ 

HIDDEN_SURF 

CLRSCRN 

2 . Subroutine INPUT 

A. Function 

Reads in or creates the initial data required for 
the elevation data file and the image data file. 

B. Input 

NAME - name of input data file 

ELFILE - elevation data file name 

EL_SIZ - elevation data number of columns 

EL_REC - elevation data number of rows 

IMFILE - image data file name 

IM_SIZ - image data number of columns 

IM_REC - image data number of rows 



ILAD, 


I LAM, I LAS - northwest 


latitude 


of 


image 


file 


ILOD, 


ILOM,ILOS - northwest 


longitude 


of 


image 


file 


FLAD, 


FLAM,FLAS - southeast 


latitude 


of 


image 


file 


FLAD, 


FLAM,FLAS - southeast 


longitude 


of 


image 


file 


TITLE 


- 50 character title placed 


on 


the 


image 



screen 

C . Output 

Same as input 

D. Calling routines 
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S0NAR3D 



E. Called routines 
CLRSCRN 

3 . Subroutine READ_ELEV 

A. Function 

Read in the elevation data file from disk 

B. Input 
None 

C . Output 

ILATD, ILATM, ILATS - southwest latitude of the 
elevation data file 

ILOND, ILONM, ILONS - southwest longitude of the 
elevation data file 

D_LAT - latitude grid size in seconds 

D_LONG - longitude grid size in seconds 

lENDM ~ number of rows of the elevation data 

lENDN - number of columns of the elevation data 

RELEV (*,*) - array containing the elevation data 

points 

D. Calling Routines 
SONAR3D 

E. Called routines 
None 
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4 . Subroutine READ_IMAGE 

A. Function 

Read in the image data file from disk 

B. Input 

IM_SIZ - number of columns of the image data 
IM_REC “ number of rows of the image data 

C . Output 

IMAGE (*,*) - array containing the image data 

D. Calling routines 
S0NAR3D 

E. Called routines 
None 

5. Subroutine TER_XYZ 

A. Function 

Converts each elevation point to it's geo- 
rectangular XYZ coordinates and calculates the 
corresponding image row/column coordinates for each 
point 

B. Input 

LATD, LATM, LATS - reference elevation latitude 
ILOND, LONM, LONS - reference elevation longitude 
D_LAT - latitude grid size in seconds 
D_LONG - longitude grid size in seconds 
RELEV (*,*) - array holding the elevation data 
lENDM - number of rows of the elevation data 



62 



lENDN - number of columns of the elevation data 



I LAD, I LAM, I LAS - reference northwest image latitude 
ILOD , ILOM , ILOS - reference northwest image 
longitude 

FLAD, FLAM, FLAS - reference southeast image latitude 
FLOD , FLOM , FLOS - reference southeast image 
longitude 

IM_SIZ - number of columns of the image data 
IM_REC - number of rows of the image data 

C. Output 

IA(*) - array containing the image column 

coordinates 

JA(*) - array containing the image row coordinates 
XYZ(*,3) - array holding the (X,Y,Z) coordinates of 
each elevation point 

D. Calling routines 
SONARS D 

E. Called routines 
C0NV2SEC 
DMS2XYZ 

6. Subroutine IM_REFAVG 

A. Function 

This routine forms the panel array containing the 
nodes of the triangles and the average gray shade. 

B. Inputs 
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IA(*) - array containing the image column 

coordinates 

JA(*) - array containing the image row coordinates 
lENDM - number of rows of the elevation data 
lENDN - number of columns of the elevation data 
IMAGE (*,*) - array containing the image data 

C. Outputs • 

PL_AR(*,5) - array containing the nodes of the 

triangle and the gray shade 

D. Calling routines 
S0NAR3D 

E. Called routines 
None 

7 . Subroutine OBS_LOC 

A. Function 

Calculates the observer location in (X,Y,Z) 
coordinates given the latitude and longitude. 

B. Input 

LATD, LATM, LATS - Latitude of observer location 
LOND, LONM, LONS - Longitude of observer location 
HEIGHT - Observer elevation in meters 
CTS - Observer course in degrees 

C . Outputs 

OBS_LOC(*) - array holding observer location 
parameters 
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X1,Y1,Z1 - Geo-rectangular coordinates of observer 
location 

X2,Y2,Z2 - Geo-rectangular coordinates of second 

observer point along LOS 

D. Calling routines 
S0NAR3D 

E. Called routines 
DMS2XYZ 

8. Subroutine M_ORIEN 

A. Function 

Determine the orientation of the image plane and 
calculate the [M] matrix parameters 

B. Inputs 

X1,Y1,Z1 - Geo-rectangular coordinates of observer 
location 

X2,Y2,Z2 - Geo-rectangular coordinates of second 

observer point along LOS 

C. Outputs 

M(3,3) - [M] matrix parameters 

D. Calling routines 
S0NAR3D 

E. Called routines 
None 
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9. Subroutine NEW_IJ 

A. Function 

Calculate the image plane coordinates of all object 
points which are in front of the field of view. 
Also mark as HIDDEN any nodes behind the image 
plane 

B. Inputs 

X1,Y1,Z1 - Geo-rectangular coordinates of the 

observer location 

X2,Y2,Z2 - Geo-rectangular coordinates of the 

second observer point along LOS 

ITOT - Total number of elevation points 

XYZ(*,3) - array holding the (X,Y,Z) coordinates of 

each elevation point 

FOCUS - Focal length 

M(3,3) - [M] matrix parameters 

C. Outputs 

IMAX(*) - X image coordinate 

IMAY(*) - y image coordinate 

DEPTH (*) - Distance from observer of each elevation 
point 

HIDDEN (*) - Flag for points hidden behind image 

plane 

D. Calling routines 
S0NAR3D 
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E. Called routines 



None 

10. Subroutine XY2IJ 

A. Function 

Convert the image (x,y) coordinates to screen (i,j) 
coordinates 

B. Inputs 

IMAX(*) - X image coordinate 

IMAY(*) - y image coordinate 

ITOT - Total number of elevation points 

XIMA_MAX - image frame size (x direction) 

YIMA_MAX - image frame size (y direction) 

HIDDEN (*) - Flag for points hidden behind image 

plane 

C. Outputs 

IA(*) - screen x coordinate 
JA(*) - screen y coordinate 

D. Calling routines 
SONARS D 

E. Called routines 
AFFIN 

11. Subroutine AFFIN 
A. Function 
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Calculates the coefficients for the AFFIN transform 
from image coordinates to screen coordinates 

B. Inputs 

XIMA_MAX - image frame size (x direction) 

YIMA_MAX - image frame size (y direction) 

C . Outputs 

A1 , A2 , Bi , B2 , Cl , C2 - parameters for the AFFIN 
transform 

D. Calling routines 
XY2IJ 

E. Called routines 
None 

12. Subroutine HIDDEN_SURF 

A. Function 

Remove hidden surfaces and plot to screen 

B. Inputs 

IA(*) - screen x coordinate 
JA(*) - screen y coordinate 

lENDM - Number of rows of the elevation data 
lENDN - Number of columns of the elevation data 
DEPTH (*) - Distance from the observer of each 

elevation point 

PL_AR(*,5) - array containing the nodes of the 

triangle and the gray shade 
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HIDDEN (*) - Flag for points hidden behind the image 
plane 

VD_ID - Virtual display identifier 
WD_ID - Window identifier 

C . Output 
None 

D. Calling routine 
SONAR3D 

E. Called routine 
QUICKSORT 
UISDC$ERASE 
UIS$SET_WRITING_INDEX 
UISDC$PLOT 

13 . Subroutine QUICKSORT 

A. Function 

Sort the panels based on maximum distance from the 
observer 

B. Input 

ARRAY(*,5) - array to sort 
KEY(*) - index to main array 
COUNT - number of elements to sort 

C . Output 

ARRAY(*,5) - original array 

KEY(*) - index to main array (revised order) 
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D. Calling routine 
HIDDEN_SURF 

E. Called routine 
None 

14 . Subroutine CLRSCRN 

A. Function 

Clear the screen 

B. Input 
None 

C . Output 
None 

D. Calling routine 
S0NAR3D 

INPUT 

E. Called routines 
None 

15. Subroutine C0NV2SEC 

A. Function 

Convert DEG, MIN, SEC coordinates to seconds 

B. Input 

DEG, MIN, SEC - Latitude or longitude in deg, min, sec 
format 

C . Output 

SEC - Total number of seconds 
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D. Calling routine 
TER_XYZ 

E. Called routine 
None 

16. Subroutine DMS2XYZ 

A. Function 

Convert DMS data to (X,Y,Z) coordinates 

B. Input 

LATD, LATM, LATS - Latitude in deg, min, sec format 
LON D, LONM, LONS - Longitude in deg, min, sec format 
HEIGHT - elevation with respect to sea level in 
meters 

C . Output 

X,Y,Z - (X,Y,Z) coordinates of input point 

D. Calling routine 
TER_XYZ 
OBS_LOC 

E. Called routine 
None 
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APPENDIX B. PROGRAM LISTING 



PROGRAM SONAR3D 

C ************************************************************ 

C 

C THIS PROGRAM OCX^STFUCTS IHE ELEVATICXI FHE AND IMAGE 

C FUE THAT IS REQUIRED BY THE 3-D PROGRAM. MAXIMUM 

C EI£VAnCXI ARRAY SIZE IS 70 RCWS X 70 OOIDMNS. MAXIMUM 

C IMAGE SIZE IS 1024 X 1024. 

C 

C ************************************************************ 
IMPLICIT INTEGER (A-Z) 

get graphics libraries 

INCLUDE 'SYS$LrBRARY:UISUSRDEF' 

INCLUDE 'SYS$LrBRARY:UISENTRY' 

INCLUDE 'SYS$LEBRARY:UISMSG' 

C 

C the follcwing variables define the elevation file inaxiinum size 
C and the image file maxiraum size. In order to simplify altering 
C of these values, they have been groijped together. Change the 
C values in the DATA statement, and in the definition statements 
C v^ch immediately follow, 
c 

GCMra /ELEV/ROW_ELEV,OOL_ELEV,TOT_EIEV,NPLANES 
OCX4MON /IMAGE/ROW_IMAGE,OOL_IMAGE 

DATA RCW_EIEV, OOL_EIEV,TOT_EIEV,NPIANES/70 , 70 , 4900 , 10000/ 
DATA RCW_IMAGE,OOL_IMAGE/1024,1024/ 

C 

BYTE IMAGE (1024, 1024) 

INTEGER HIDDEN (4900) ,IA(4900) ,JA(4900) 

INTEGER PL_AR(10000,5),DEPIH_KEY(10000) 

REAL REIEV(70,70) 

REAL*8 DEPIH(4900) ,XYZ(4900,3) ,IMAX(4900) ,IMAY(4900) 

C 

CHARACTER EIFILE*13,IMFII£*13,FILE_NAME*13,ANS*l,TrTIE*50 
C 

INTEGER IENm,IENIM,ITOT,IIATD,IIAIM,IIATS 

INTEGER II£M),IL:»^M,II0NS 

INTEGER D_IAT,D_IONG,IATS,IONS,OBSIOC(8) 

INTEGER EL_SIZ,EL_REC,IM_SIZ,IM_REC 
INTEGER IIAD,IIAM,HAS,IIOD,IIOM,IIOS 
INTEGER FlJ^,FIAM,FIAS,FL3D,FICM,FIOS 
C 

REAL*4 I_VECT0R 

REAL*8 X1,Y1,Z1,X2,Y2,Z2,M(3,3) 
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REAL*8 FOCUS, XFRAME,YFRAME 
PEAL*8 HEIGHT, FWR, SCALE 
C 

C set i:?) graphics environment 
C 

C create a color map with 251 entries 
C 

DATA VCM_SIZE/251/ 

V(ll_II><nS$CREATE_OOIDR_MAP (VCM_SIZE) 

C 

C make a virtual display 16.8 cm x 16.8 cm 
C 

VD_IIMJIS$CE^EATE_DISPIAY (0 . 0 , 0 . 0 , 512 . 0 , 512 . 0 , 16 . 8 , 16 . 8 , VCM_ID) 
C 

C fill in the color map 
C 

DO 50 1=0,250 
I_VECrOR=I/250. 

CALL UIS$SET_INrENSITY (VD_ID, I , I_VECIOR) 

50 CONTINUE 

C 

C set the background color and fill pattern 
C 

CALL UIS$SET_INrENSr]:Y (VD_ID, 0,0.0) 

CALL UIS$SET_FCWr(VD_ID, 1,1, 'UIS$FILL_PATrERNS ' ) 

CALL UIS$SET_FILL_PATIERN(VD_ID,l,l,PArr$C_POREGROUND) 

C 

C start main routine, read irput 
C 

CALL INIUT(ELFH£,EL_SIZ,EL_REC, 

IMFIIE, IM_SIZ , IM_REC/ HAD, IIAM, IIAS, IIDD, Iim, IIDS , 

FIAD, FIAM, FIAS , FIDD, Fim, FIDS , TITLE) 

C 

EL_REC=EL_REC+1 

IM_SIZ=IM_SIZ/4 

C 

C open the input data files 
C 

OPEN(UNIT=l,FII£=EIEH£,STAIUS='OID' ,AOCESS=' DIRECT' , 
REC0RDSIZE=EL_SIZ , MAXREO=EL_REC) 

OPEN (UNIT==4,FIIE=IMFII£,STAIUS=' OLD' ,AOCESS=' DIRECT' , 
REO0RDSIZE=IM_SIZ , MAXREO=IM_REC) 

C 

CALL READ_ELEV(IIATD,ILAIM,IIATS,IICND,ILa^M,ILa^S, 

D_IAT, p_D:»rc, REIEV, lENEN, lENEM) 

CALL READ_IMAGE ( IMAGE , IM_SIZ , IM_REC) 

CALL TERJCYZ(IIATD,ILAIM,IIATS,II£M),IIDNM,Iia^S, 

D_IAT, D_L:M5, REIEV, lENEM, lENEN, XYZ , lA, JA, 

IIAD, HAM, IIAS , IICD, IICM, IIDS , 

FIAD, FIAM, FIAS , FIDD, FILM, FLOS , IM_SIZ , IM_REC) 

C 

CALL IM_REFAVG(IA,JA, IMAGE, IENEM,IENEN,PL_AR) 
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YFRAME=29000.0 

Xn?AME>=35000.0 

FIIE_NAME=' ' 

SAVE_FIAG=0 

P0CUS=.050 

H^l.O 



SCAIE=1.2 



nOT=IENCN*IENEM 

WD_II>^JIS$CKEATEj™XW (VD_ID, ' SYS$V«3RKSTATI<^ ' , TTIIE , 
0.0,0.0,512.0,512.0,16.8,16.8) 

C 

C 

60 (2ALL 0BS_IDC(X1,Y1,Z1,X2,Y2,Z2,0BSIDC) 

(2ALL M_0RIEN(M,X1,Y1,Z1,X2,Y2,Z2) 

CMi NEW_U (XI, Y1,Z1,X2,Y2,Z2, nor, XYZ, FOCUS, M,IMAX, 
IMAY , HIDDEN , DEPIH) 

70 C3UX XY2IJ(IA,JA,IMAX,IMAY,XFRAME,YFE^AME,IT0T, HIDDEN) 

CALL HIDDEN_SURF(IA,JA,IENEM,IENEN, DEPIH, PL_AR,VD_ID, 

. HIDDEN , DEPIH_KEY , WD_ID) 

C 



C type the c±)server location and menu 



C 

80 



85 

90 

91 



C 



CALL dPSCKN 

WRITE(6,*) 'OBSERVER LDCATKa^' 

WRITE(6,*) 

WRITE(6,85) 'lAT: ' ,OBSIOC(l) ,OBSIDC(2) ,OBSIOC(3) 
WRITE(6,85) 'lOTC: ' ,OBSIOC(4) ,OBSIDC(5) ,OBSIOC(6) 
WRITE(6, 90) 'HEIGHT; ' ,OBSIDC(7) , ' HEADING: ',OBSIDC(8) 
WRITE(6,91) 'MAGNIFICATIC»I: ',FWR, ' YSCAIE: ', SCALE 
FORMAT (A8, 14, 13, 13) 

FORMAT (A9 , 17 , A12 , 13 ) 

FORMAT (A16 , F4 . 1 , A12 , F4 . 1) 

WRITE(6,*) 

WRITE(6,*) 'SELECT OF THE FOLLCWING' 

WRITE (6,*) ' (M)agnification' 

WRTTE(6,*)' (Y)scale factor' 

WRITE (6,*)' (O)bserver Location' 

WRTTE(6,*)' (S)ave linage' 

WRTTE(6,*)' (E)xit' 



READ (5, 100) ANS 
100 FORMAT (Al) 

C 



IF (ANS .NE. 'M' .AND. ANS .NE. 'm') GO TO 105 
WRTTE(6,*) 'ENTER MAGNIFICATION FACTOR' 
READ(5,*)IWR 

102 XFRAME=35000 . 0/IWR 

YFRAME=35000. 0/ (IWR*SCALE) 

GO TO 70 



C 

105 IF (ANS .NE. 'Y'.AND. ANS .NE. 'y') GO TO 110 
WRTTE(6,*) 'ENTER Y-SCAIE FACTOR' 
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PEAD(5,*) SCALE 
GO TO 102 
C 

no IF (ANS .NE. 'O'. AND. ANS .NE. 'o') GOTO 115 
GO TO 60 
C 

115 IF (ANS .NE. 'S' .AND. ANS .NE. 's') GO TO 118 

IF (FH£_NAME .EQ. ' ') HiEN 

WRITE(6,*) 'ENTER IHE NAME OF IHE IMAGE FILE' 
READ (5, 116) FH£_NAME 
ENDIF 

SAVE_FIAG==SAVE_FIAG+1 

CALL IMAGE_SAVE(MD_ID,FIIE_NAME,SAVE_FIAG) 

GO TO 80 

116 FORMAT (A13) 

118 IF (ANS .NE. 'E'.AND. ANS .NE. 'e') GO TO 80 
C 

120 IF (SAVE_FLAG .GT. 0) THEN 

WRITE (11) SAVE_FLAG 
ENDIF 
CLOSE (1) 

CLOSE (4) 

CLOSE ( 10 , STATUS= ' SAVE ' ) 

CLOSE ( 11 , STATUS= ' SAVE ' ) 

CALL UIS$DELETE_DISPLAY(VD_ID) 

CALL CLRSCRN 
END 



75 



I 



c 

C ****icirkirkirk*ic**************************icirkitirkirk**ic***** 

C 

SUBROUTINE INH7r(EiraE,EL_SIZ,ELJ?EC,IMFIIE, 

IM_SIZ , IM_REC, IIAD, IIAM, IIAS , HDD, II£M, nos , 

FLAD, FIA5 , FIOD, FICM, FIOS , TITLE) 

C 

C INR7IS = 

C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 

c ounuis = 
c 

C **************: 

c 

IMPLICIT 
C 

CHARACTER ELETIE*13,IMFII£*13,TITIE*50,ANS*1,NAME*13 
C 

INTEGER EL_SIZ,EL_REC,IM_SIZ,IM_REC 
INTEGER ILAD, HAM, HAS, FLAD, FLAM, FIAS 
INTEGER HDD, IICM,HOS,FIOD, FLOW, FLOS 
C 

C start routine here 
C 

5 CALL CLRSCRN 

WRITE(6,*) 'Do you wish to create an iiput data file (Y/N)? ' 
READ (5, 100) ANS 
100 FORMAT (Al) 

IF (ANS .EQ. 'Y' .OR. ANS .EQ. 'y') GO TO 200 
IF (ANS .EQ. 'N' .OR. ANS .EQ. 'n') GO TO 500 
GO TO 5 
C 

C create irput file 
C 

200 CALL CLRSCRN 

WRITE(6,*) 'Iiput the filename for the input file(*.dat) : ' 
READ(5,205) NAME 

OPEN (UNIT=1,NAMB=NAME,STAIUS='NEW', 

ACCESS= ' SEQUENTIAL' , FORM= ' FORMATTED ' ) 

C 

C get elevation file data 
C 



ELFHE elevation data file name 
EL_SIZ # column of ELFHE 

EL_REC # rows of ELFHE 

IMFHE image data file name 
IM_SIZ # columns of IMFHE 

IM_REC # rows of IMFHE 

HAD, HAM, HAS 
HOD,HCM,HOS 
FLAD,FIAM,FLAS 

FIOD,FLCM,FLOS Southeast comer of image file 
NAME file name for the input data file 
TITLE header name for the image 



Northwest comer of image file 



NO® 



It ************:!:**************************** 



INTEGER (A-Z) 
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WRITE(6,*) 'Enter the elevation data file name: ' 
READ(5,205) EIJTLE 

WRITE(6,*) 'Elevation file number of columns: ' 

READ(5,*) EL_SIZ 

WRITE(6,*) 'Elevation file number of rows: ' 

READ(5,*) EL_REC 

get image file data 

WRITE (6,*) 'Enter the image data file name: ' 

READ(5,205) IMFTIE 

WRITE(6,*) 'Image file number of columns: ' 

READ(5,*) IM_SIZ 

WRITE(6,*) 'Image file number of rows: ' 

READ(5,*) IM_REC 

northwest comer of image latitude and longitude 

WRITE (6,*) 'Enter the North-West lat/long of image' 
WRITE(6,*) 'Irpjt the latitude (DEG,MIN,SEC) : ' 

READ (5 , * ) IIAD, IIAM, IIA5 

WRITE(6,*) 'Infxit the longitude (DEG,MIN,SEC) : ' 

READ(5 , * ) HDD, H£M, ILOS 

southeast comer of image latitude and longitude 

WRITE (6,*) 'Enter the South-Ecist lat/long of image' 
WRITE(6,*) 'Irput the latitude (DEG,MIN,SEC) : ' 

READ(5,*) P1AD,FLAM,FIAS 

WRITE(6,*) 'Irput the longitude (DEG,MIN,SEC) : ' 

READ(5,*) FIDD,FIDM,FIJDS 

WRITE (6,*) 'Irpjt a title for the image (50 char max.) : ' 
READ(5,206) TITLE 

write data to data file 



C 

205 

206 
207 
209 



WRITE(1,205) 
WRITE(1,207) 
WRITE(1,207) 
WRITE (1,205) 
WRITE(1,207) 
WRITE(1,207) 
WRITE(1,209) 
WRITE(1,210) 
WRITE(1,209) 
WRITE(1,210) 
WRITE(1,206) 



ELFTIE 

EL_SIZ 

EL_REC 

IMFTIE 

IM_SIZ 

IM_REC 

HAD, HAM, HAS 
HDD,n£M,HDS 
FIAD,F1AM,F1AS 
F10D,F1£M,FIDS 

TTTTF 



PQPMAT(A13) 

F0iPMAT(A50) 

P0RMAT(I3) 

FOFMAT(3I3) 
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210 FOEMAT(I4,2I3) 
GO TO 600 



read an input file 



500 



C 

600 



CALL dl^SCRN 

VZRITE(6,*) 'Enter the name of the input data file; 
PEAD(5,205) NAME 

OEEN(lirrr=l,NAMB^^AME,STAIUS='OID' , 

AOCESS= ' SEQUENITJUj* , POBJ^ • POEMATTED ' ) 



PEAD(1,205) 

PEAD(1,207) 

PEAD(1,207) 

PEAD(1,205) 

PEAD(1,207) 

READ(1,207) 

READ(1,209) 

READ(1,210) 

READ(1,209) 

READ(1,210) 

PEAD(1,206) 



FTFTTF! 

EL_SIZ 

EL_REC 

IMFIIE 

IM_SIZ 

IM_EEC 

HAD,IIAM,IIAS 

HDD,I]XM,IIOS 

F1AD,FLAM,FIAS 

FIDD,FIXM,FIDS 

TITIE 



I 



CLOSE (UNTT=1) 

PEIURN 

END 
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4 : ^ 



c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



SUERCUITNE READ_EIEV(I]jm),IIAIM,IIAIS,IKWD,IlJ:X^,IlJ:»^S, 
D_IAT, D_I£NG, PEIEV, lENDN, lENEM) 



THIS SUEROUTTNE READS THE ELEVATION FIIE AND PLACES THE DATA 
INTO A REAL 50X50 ARRAY RELEV() . 

INH7IS = NONE 



OUTEUT = HATD 


southwest 


IIATM 


southwest 


HATS 


southwest 


HOLD 


southwest 


nCML 


southwest 


HONS 


southwest 


D LAT 


grid size 


D icm 


grid size 


RELEVO 


elev data 


lENEM 


# rows of 


lENEN 


# columns 



latitude of elev data (degrees) 
latitude of elev data (minutes) 
latitude of elev data (seconds) 
longitude of elev data (degrees) 
longitude of elev data (minutes) 
longitude of elev data (seconds) 
in minutes of latitude 
in minutes of longitude 

elev data 
of elev data 



******************************************************************* 



IMPLLCIT INTEGER (A-Z) 

C 

OCMMDN /ELEV/R0W_EI£V,00L_E1£V 
INTEGER IENEN,IENEM,D_1AT,D_IJCM3 
INTEGER IIATD,IIATM,I]ATS,IL»ro,IlJDNM,IICNS 
C 



REAL RELEV (R0W_ELEV, 00L_ELEV) 

C 

READ(1,REC=1) ILATD, IIATM,IIATS,inMD,mMl, IimS, 
D_IAT, D_IiDNG, lENDN, lENEM 
C 

DO 10 M=1,IENEM 

READ(1,RE0^1) (RELEV (M,N) ,N=1,IENEN) 

10 CONTINUE 

RETURN 
END 
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************************************************************ 

SUBROUTINE READ_IMAGE (IMAGE, IM_SIZ,IM_FEC) 

C 

C THIS SUBROUTINE READS IHE IMAGE DATA INTO AN ARRAY. 

C INIUTS = IMJ5IZ number of columns of image 

C IM_REC number of news of image 

C 

C OUnUT = IMAGE 0 image gray level file 

C 

Q'kiciciciciciciciciciciciciticicicicicicicicicicicicicicicicicicicicicicicicicicicicicicic'k'k'kicicicicicicicicicicicic 

IMPLICIT INTEGER (A-Z) 

OCMWa^ /IMAGE/ROW_IMAGE,OOL_IMAGE 
C 

BYTE IMAGE (ROW_IMAGE,OOL_IMAGE) 

C 

INTEGER IM_SIZ,IM_REC 
C 

DO 10 IR=1,IM_REC 

READ(4,REC=IR) (IMAGE(IR,IC) ,IC=1,IM_SIZ*4) 

10 OC»ITINUE 
RETURN 
END 
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*************************************************************** 



SUBROOTINE TER_XYZ(UOT),1J^,IAIS,IIj:M),IDNM,DDN 
. RELEV, IE^O^,IE^CN,XYZ,IA,JA, 

. IIAD,IIAM,IIAS,IIDD,ILCM,IIDS, 

. FIAD,FU^,FIJ^,FIDD,Fim,FIJDS,IMjSIZ,IM_REC) 

THIS SUBROUTINE CONVERTS EACH ELEVATK^ POINT TO ITS EMS 
BQUIVAIENT. IT USES THE FACT THAT EACH POINT REPRESENTS 
AN EQUAL CHANGE IRCM THE IA5T. THE REFERENCE IAT/KM3 
AND THE DELTA lAT/ILMS ARE PASSED IN THE SUBROUTINE CALL. 

. . .ASSUMES COLOR CCMPLETEL^ OVERLAPS THE EIEV DATA 

. . .NEED TO KNOW THE lAT/IONG OF THE OOim IMAGE 

. . .Also CALCS THE IA() , JA() COORDINATES FOR EACH EIEV FT 

EIEVATICH DATA 

INPUTS lATD - REFERENCE lAITTUDE (DEGREES) 
lAIM - REFERENCE lAITTUDE (MINUTES) 
lATS - REFERENCE lAITTUDE (SECONDS) 

IL»ID - REFERENCE LONGITUDE (DEGREES) 
liML - REFERENCE KMGITUDE (MINUTES) 

I£»LS - REFERENCE LONGITUDE (SECONDS) 

D_IAT - DISTANCE BETWEEN ROWS 
D_iarc - DISTANCE BETWEEN COLUMNS 
REIZVO - ELEVATION DATA FILE 
lENEM - NUMBER OF RCWS 
lENEN - NUMBER OF COLUMNS 

OUTEUTS XYZ - OUTPUT DATA FILE 

******************************************************************* 
IMPLLCrT INTEGER (A-Z) 

CCMMON /ELEV/ RCW_EIEV,OOL_EIEV,TOT_EIEV 

INTEGER IA(TOT_EIEV) , JA(TOT_EIEV) 

INTEGER IATD,IAIM,IAIS,ILOND,IL:Ml,D_IAT,D_ia^ 

INTEGER LOND,U:ML,imS 
INTEGER IIAD,HAM,IIAS,ILOD,IICM,ILOS 
INTEGER FIAD, FLAM, FLAS,FLOD,FLCM, FLOS 
INTEGER IM_SIZ,IM_REC 

REAL RELEV (RCW_EIEV,OOL_EIEV) 

REAL*8 HEIGHT, X,Y,Z,XVrZ(TOT_EIEV, 3) 

convert to seconds 

CALL OC»W2SEC(IIAD,IIAM,IIAS) 

CALL CC»^V2SEC(IL0D,II£M,IL0S) 

CALL OCXW2SEC (FIAD, FLAM, FLAS) 

CALL CC»^V2SEC(FL0D,F1CM,FL0S) 
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DELLAI^IIAS-FLAS 
DELLDN=FIDS-IK)S 

set flag for proper hemisphere 

IF (HAD .GE. 0) UiEN 
IAT_FLAG =1 
ELSE 

IAT_FIAG=-1 
ENDIF 

IF (HOD .GE. 0) THEN 
ION_FEAG=i 
ELSE 

IDN_FIAG =-1 
ENDIF 

find initial seconds of elev. data 

CALL OC»^V2SEC(IATD,IAIM,IATS) 

CALL O»JV2SEC(imC),I0NM,Ii:»JS) 

set 1 : 5 ) for 1 st increment 

IAIS=IATS-D_IAT 

hons=ions-d_i£m; 

start routine here 
IR=0 

DO 10 M=1,IENEM 
IAIS=IAIS+D_IAT 
IONS=in»IS 
DO 20 N=1,IENEN 
IR=IR +1 

IONS=rONS+D_IONG 
HEIGHI^RELEV(M,N) 

R0W=INr ( (IM_REC-1) * (IIAS-IATS)/DELLAT) +1 
OOL=INr ( ( ( IM_SIZ*4 ) -1) * (K^JS-IIOS) /DELLIMT) +1 

the row and column coordinates are 1-512 

IA(IR)=<X)L 
JA(IR)=RCW 

CALL EMS2XYZ ( 0 , 0 ,IAIS, 0 , 0 ,lJ:»^,HEIGHr,X,Y,Z) 
XYZ(IR,1)=X 
XYZ(IR,2)=Y 
XYZ(IR,3)=Z 
20 CONTINUE 

10 CONTINUE 
RETURN 
END 
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******************************************************************** 

SUBROUTINE CMS2XYZ (IATD,IAIM,IATS,I0ND,n»^M,Ii»TS,HEIC3n’,X, Y,Z) 
C 



C IHIS SUBROUTINE CONVERTS EMS DATA TO X,Y,Z. 

C 



c 


INIUTS = lAID 


latitude in degrees 




c 


lAIM 


latitude in minutes 




c 


LATS 


latitude in seocsnds 




c 


lOND 


longitude in degrees 




c 




longitude in minutes 




c 


IONS 


longitude in seconds 




c 

c 


HEIGHT 


hei^t above sea level 


(meters) 


c 


OUTFUT = X,Y,Z 


XYZ coordinates of EMS 


point 



C 

c ********************************************************************* 

IMPLrCTT INTEGER (A-Z) 

C 

INTEGER IATO,IAIM,IATS,IOND,D»W,I£»iS 
C 

REAL*8 IHI,IAMnA,N,X,Y,Z, HEIGHT 
REAL*8 PI,C1,C2,C3,E_SQUARE,A,RADIAN 
C 

PARAMETER (PI=3 . 14159265358793238) 

PARAMETER (Cl=180 . , C2=60 . , C3=3600 . ) 

PARAMETER (E_SQUARE=0 . 006768658 ,A=6378206 . 4 ) 

C 

RADIAI^PI/Cl 

C 

IF (lATD .GE. 0) THEN 
FHI= (lATDf IADVC2+IATS/C3 ) *RADIAN 
ELSE 

FHI= (IATD-IADVC2-IATS/C3 ) *RADIAN 
END IF 
C 

IF (mro .GE. 0) THEN 
IAMDA= (lONDf ION^VC2+IONS/C3 ) *RADIAN 
ELSE 

IAMDA= (IOND-lJOi^C2-IONS/C3 ) *RADIAN 
END IF 
C 

N=^A/SQRT(1-E_SQUARE*SIN(IHI) *SIN(FHI) ) 

X=(N+HEIGHT) *OOS (FHI) *00S (lAMDA) 

Y=(NHiEIGHT) *OOS (FHI) *SIN(IAMDA) 

Z= (N* (1-E_SQUARE) +HEIGHT) *SIN (FHI) 

RETURN 

END 
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c 



c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



SUEROUITNE IMJ?EEAVG(IA,JA, IMAGE, IEtOd,IEtrai,PLJ^) 



THIS SUBROUTINE CXX^SIRQCIS THE GBCMEERY FIIZ 
THAT CONTAINS THE THREE NODAL POINTS THAT MAKE UP 
AN IMAGE PLANE AND THE GREY SCAIE VALUE ASSIGNED TO 
THAT PLANE. 



INPUTS = IMAGEO 
lENEM 
lENEN 
IA() 
JA() 



image gray level file 

# rcws of elev ciata 

# cx5lurtins of elev data 
image oolimm index 
image rcM index 



OUnUT= PL_AR() cirray holding nodes and 

gray value for each panel 



*********************************************************** 



IMPLICIT INTEGER (A-Z) 

C 

CCMM:»T /EIEV/RCW_EIE7,GOL_EIEV,TOT_EIEV,NPIANES 

ccmm:»t /image/row_image,gol_image 

c 



BYTE IMAGE (RCW_IMAGE,COL_IMAGE) 
C 



INTEGER IA(TOT_EIEV) ,JA(TOT_EIEV) ,PL_AR(NPIANES,5) 
INTEGER IY,M,N,IGREY1,IGREY2,L,L1 
INTEGER IENEM,IENEN,NODE_A,NODE_B,NODE_C,NODE_D,IR 
INTEGER IGOUNT1,ICOUNT2,ITOT1,ITOT2 



REAL*8 SLOPE, YINT 
C 

DO 90 H^1,IENEM-1 
IL=(IR-1) *IENEN 
DO 80 N=1+IL,IENCN-1+IL 
NODE_A=N 
NODE_B=Nfl 
NODE_C=NfIENEN 
NODEJMraDE_C+l 

SIOPE= ( JA (NODE_C) - JA (NODE_B) ) *1 . 0/ (lA (NODE_C) -lA (NODE_B) ) *1 . 0 
YINT=1 . 0*JA(NODE_B) -SIOPE*IA(NODE_B) 

C 

C new knew the slope and y-intec^jt of the line forming the triangle 
C 



non=o 

noT2=o 

IC0UNT1=0 

IC0UNT2=0 

DO 70 M=IA(NODE_B),IA(NODE_A),-l 
IY=SIOPE*MfYINT 
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DO 50 1XIA(N0DE_B) 

IF (IMAGE(L,M) .ET. 0) THEN 
IT0T1=ITOT1+IMAGE (L, M) +256 
ELSE 

IT0T1=ITC)T1+IMAGE (L,M) 

END IF 

ioouNri=iooui7ri+i 

50 OXOTHUE 

IF (M. IT . lA (NODE_B) ) IHEN 
DO 60 Lr=IY,JA(NODE_C) 

IF (IMAGE(L,M) .IT. 0) IHEN 
IT0T2=riOT2+IMAGE (L,M) +256 
EISE 

IT0T2=riOT2+IMAGE (L, M) 

END IF 

IOOUNr2=IOOUNT2+l 
60 OCWTINUE 

END IF 

70 OCmiNUE 

L1=(N-(IR-1))*2 
IGKEYl=ITOTl/IOOUNri 
IGREY 2 =riOr 2 /IOOUNT 2 
PL_AR (Ll-1 , 1) =N0DE_A 
PL_AR (Ll-1 , 2 ) =N0DE_B 
PL_AR(L1-1, 3)=N0DE_C 
PL_AR(L1-1, 4) =IGREY1 
PL_AR(L1 , 1) =f«)DE_B 
PL_AR (LI , 2 ) =NODE_D 
PL_AR (LI , 3 ) =NODE_C 
PL_AR(L1, 4) =IGREY2 
80 OCmiNUE 
90 OCmiNUE 
RETURN 
END 



85 



VOOO HO oooo ooooooooooooooo 



************************************************************ 
SUBROUTINE OBS_IDC(X1,Y1,Z1,X2,Y2,Z2,OBSIOC) 

ims SUBROUTINE CALOJIATES TOE NEW OBSERVER XYZ 
lOCATiai FROM DESIRED lAT. AND I£»JG. INIUTS. 

INHJTS = NONE 

OUTEUTS = X1,Y1,Z1 observer location 

X2,Y2,Z2 observer locaticxi direction 

OBSIOCO observer location array 

************************************************************ 
IMPLICrT INTEGER (A-Z) 

CHARACTER ANS*1 

INTEGER LATD,IATO,IATO,im),I£m,ia^S,CTS 
INTEGER OBSIOC(8) 

REAL*8 X1,Y1,Z1,X2,Y2,Z2,HEIGHT 

read in last values of observer location 

IATD=OBSIOC(1) 

LAmOBSIOC(2) 

IATS=OBSIOC(3) 

IOO=OBSIOC(4) 

IJDNM=OBSIOC(5) 

imS=OBSIOC(6) 

HEIGHI^BSIOC(7) 

dS=OBSIOC(8) 

F1AG=0 

CALL CIRSCRN 

WRITE (6,*) 'OBSEEA/ER DDCATIQN' 

WRTTE(6,*) 

WRTTE(6,8) 'lAT: ' ,OBSIOC(l) ,OBSIOC(2) ,OBSIOC(3) 
WRTTE(6,8) 'LCX^: ' ,OBSIOC(4) ,OBSIOC(5) ,OBSIOC(6) 

WRITE (6, 9) 'HEIGHT: ' ,OBSIOC(7) , ' HEADING: ',OBSIOC(8) 
FORMAT (A8, 14, 13, 13) 

FORMAT (A9 , 17 , A12 , 13 ) 

WRTTE(6,*) 

IF (FLAG .EQ. 1) THEN 
WRTTE(6,*) 'ARE TOERE ANY OORRECnO^S ' 

WRTTE(6,*) 

ENDIF 

WRTTE(6,*) 'SELECT GtJE OF TOE FOLLOWING' 

WRTTE(6,*)' (1) Latitude' 

WRTTE(6,*)' (2) Longitude' 

WRTTE(6,*)' (3) Hei^t' 
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WRITE(6,*)' (4) Heading* 

WRITE(6,*)' (5) All* 

WRITE(6,*)' (6) None* 

C 

READ(5,100) ANS 
100 FORMAT (Al) 

IF (ANS .EQ. *6') GO TO 160 
C 

IF (ANS .NE. *1* .AND. ANS .NE. '5') GO TO 60 
WRITE(6,*) 'INIUr OBSEE^VER lATTIUDE -DEGREES: ' 

READ(5, *) lAID 

WRITE(6,*)' -MINUTES: ' 

READ(5,*)IAIM 

WRITE(6,*)' -SECONDS: ' 

READ(5,*)IATS 
C 

60 IF (ANS .NE. '2*.AND. ANS .NE. *5') GO TO 70 

WRITE(6,*) 'INHTT OBSERVER lO^GITUDE -DEGREES: ' 

READ(5,*)I£M) 

WRITE(6,*)' -MINUTES: ' 

READ(5,*)LCM^ 

WRrTE(6,*)' -SECONDS: ' 

READ(5,*)IiONS 
C 

70 IF (ANS .NE. *3 '.AND. ANS .NE. *5') GO TO 80 

WRrTE(6,*) 'INIUT OBSEF!VER HEIGHT ^4ETERS: ' 
READ(5,*) HEIGHT 
C 

80 IF (ANS .NE. *4* .AND. ANS .NE. '5') GO TO 105 

WRTTE(6,*) 'INIUT IHE COURSE DEGREES: ' 

READ(5,*)CTS 

105 IF (ANS .IT. *1* .OR. ANS .GT. *6') GO TO 1 

C 

C load t±ie values into the obs_loc array 
C 

150 OBSIDC(l)=IATD 
OBSIDC(2)=IAIM 
OBSIDC(3)=IAIS 
OBSIJOC(4)=IJOND 
OBSIJOC(5)=L»IM 
OBSDOC(6)=IDNS 
OBSIDC(7)=4ffiIGHT 
OBSDOC(8)=CIS 
FIAG=1 
GO TO 1 
C 

C convert to XYZ coordinates 
C 

160 CALL CMS2XYZ(IATD,IAIM,IATS,mro,lJ0NM,Ii0NS,HEIGHT,Xl,Yl,Zl) 
C 

C calc new position based on course 
C 
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IF (lATD .GT. 0) UffiN 

IATS=5*00S (2*3 . 14159*CIS/360)+IATS 
ELSE 

IArS=IATS-5*SIN (2*3 . 14159*CTS/360) 

END IF 

IF (UM) .GT. 0) THEN 

IJDNS=5*SIN (2*3 . 14159*CTS/360) +KXIS 
ELSE 

L£»JS=LmS-5*SIN (2*3 . 14159*CTS/360) 

END IF 

convert new positiai to XYZ coordinates 

CALL EMS2XYZ (LAID, lAIM, LAIS, HEIGHT, X2,Y2,Z2) 

REIURN 

END 
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Q ************************************************************ 

SUBRCUITNE M_0RIEN(M,X1,Y1,Z1,X2,Y2,Z2) 

C 

C IKES RDOTINE DETEFMINES THE ANGIES OF ROTATION 

C OF THE IMAGE PLANE WITH RESPECT TO THE WORIE COORDINATE AXIES. 

C AND C2UX3JIATES THE NEW M MATRIX FOR THE NEW VIEWER LOGATTON 

C 

C INIurs = X1,Y1,Z1 observer location 

C X2,Y2,Z2 observer direction point 

C 

C OOTEUTS = M() 3D - 2D xform matrix 

C 

C *************************************************************** 
IMPUCTT INTEGER (A-Z) 

C 

REAL*8 MAGN_X,MAGN_Y,MAGN_Z 

REAL*8 X_VECX,X_VECY,X_VECZ,YJTX3C,Y_VECY,Y_VECZ 
REAL*8 Z_VECX,Z_VECY,Z_VECZ,X1,Y1,Z1,X2,Y2,Z2 
REAL*8 M(3,3) 

C 

C get the two OBS IOC points 
C and set 15) vectors 
C 

Y_VECX=X1 

Y_VECY=Y1 

Y_VECZ=Z1 

C 

C select another point along the track 
C 

Z_VECX=X1-X2 

Z_VECY=Y1-Y2 

Z_VECZ=Z1-Z2 

C 

C use the cross product of Y CROSS Z to obtain the X vector 
C 

X_VECX=( (Y_VECY*Z_VECZ) -(Y_VECZ*Z_VECY) ) 

X_VECY=( (Y_VECZ*Z_VECX) -(Y_VECX*Z_VECZ) ) 

X_VECZ=( (Y_VECX*ZJTX::Y) -(Y_VECY*Z_VB2X) ) 

MAGN_Z=SQRT( (Z_VECX**2) + (Z_VECY**2) + (Z_VECZ**2) ) 

MAGN_X=SQRT( (X_VECX**2) + (X_VECY**2) + (X_VECZ**2) ) 

MAGN_Y=SQRT ( (Y_VECX**2 ) + (Y_VECY**2) + (Y_VECZ**2 ) ) 

M ( 1 , 1) =xj/ecvm?^_x 
M ( 1 , 2 ) =X_VECY/MAGN_X 
M(l, 3) =X_VECZ/MAGN_X 
M(2 , 1) =Y_VECVMAGN_Y 
M(2 , 2) =Y_VECY/MAGN_Y 
M(2 , 3) =Y_VECZ/MAGN_Y 
M ( 3 , 1) =Z_VECVMAGN_Z 
M (3 , 2 ) =Z_VECY/MAGN_Z 
M(3 , 3) =Z_VECZ/MAGN_Z 
RETURN 
END 
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**************************************************************** 

SUBROUTINE NEW_U(X1,Y1,Z1,X2,Y2,Z2,IT0T, 

XYZ , FOCUS IMAX, 11^, HIDDEN, EffiPIH) 



IHIS SUBROUTINE OCMFUTES THE NEW lA AND JA SCREEN 
COORDINATES FROM THE GIVEN OBSERVER lOCATION. 



INPUTS 



OUTFUIS 



X1,Y1,Z1 


observer location 


X2,Y2,Z2 


observer direction point 


rroT 


# elev points 


M() 


3d - 2d xform matrix 


FOCUS 


focal length 


XYZ() 


dms coordinates of elev data 


IMAXO 


X image coordinate 


IMAYO 


y image coordinate 


HIDDEN 0 


flag for points hidden behind FOV 


DEPIHO 


dqpth of each point 



*************************************************************** 
IMPLICIT INTEGER (A-Z) 



OOMMJN /ELE^/ ROW_EIEV,OOL_EIE7,TOT_EIEV 

INTEGER rrOT,IR, HIDDEN (TOT_EI£V) 

REAL*8 X1,Y1,Z1,F0CUS,M(3,3) 

REAL*8 X,Y,Z,XIMA,YIMA,DENCM,DEPIH(TOT_EIEV) 

REAL*8 XYZ (TOT_EI£V, 3 ) , IMAX (T0T_EIEV) , IMAY (TOT_ELEV) 
REAL*8 X2,Y2,Z2,OBS_VECX,OBS_VECY,OBS_VECZ 
REAL*8 OBJ_VECX,OBJ_VECY,OBJ_VECZ,MAG_OBS,MAG_OBJ 
REAL*8 OOS_THETA,VAIUE 



VAIUE=90.0 

VAIUE=COS (VAIUE/57 . 29577951) 



DO 20 IRF=l,nOT 
X=XYZ (IR, 1) 
Y=XYZ(IR,2) 
Z=XYZ(IR,3) 

calc the c±s vec 



OBS_VECX=X2-Xl 

OBS_VECY=Y2-Yl 

OBS_VECZ=Z2-Zl 

MAG_OBS=SQRT (OBS_VECX**2-K)BS_VECY**2+OBS_VECZ**2 ) 
calc the obj_vec 
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0BJ_VECX=X-X2 
0BJ_VECY=Y-Y2 
0BJ_VECZ=Z-Z2 

MAG_OBJ=SQRT (0BJJ/ECX**24OBJ_VECY**24OBJ_VECZ**2 ) 

calc cos theta, vhere theta is angle between line of si^t 
and object point. 

<XSJIHEI»fOBS_VECX*0BJ_VECX+OBS_YEX:Y<OBJJVECY4OBSJ/ECZ*0BJ_V^ 
OOSjmETA=<X)S_THErA/ (MAG_OBS*MAG_OBJ) 

new check if (DOS (theta) > 0 => angle < 90 degrees 

IF ((DOSJEHETA .GT. VAIUE) IHEN 

DENCM=M(3 , 1) * (X-Xl) 4M(3 , 2) * (Y-Yl) 4M(3 , 3) * (Z-Zl) 

XIMA=-P0aJS* (M ( 1 , 1) * (X-Xl) 4M (1 , 2 ) * (Y-Yl) 4M ( 1 , 3 ) * (Z-Zl) ) /DENOM 
YIMA^FOCUS* (M (2 , 1) * (X-Xl) 4M (2 , 2) * (Y-Yl) 4M(2 , 3 )* (Z-Zl) ) /DENOM 
XIMA=XIMA* 1000000. 0 
YIMA=YIMA*1000000. 0 

save xima, yima in IMAX() , IMAY() array terrporarily 

IMAX(IR)=XIMA 
IMAY(IR)=YIMA 

DEPIH(rR) =SQRT( ( (Xl-X) **2) + ( (Yl-Y) **2) + ( (Zl-Z) **2) ) 

HIDDEN (IR)=0 
ELSE 

HIDDEN (IR)=1 
END IF 
20 crwriMJE 
RETURN 
END 
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***************************************************************** 



SUBROUTINE XY2U(IA,JA, 

IMAX, IMAY,XroAIffi, IT0T,Hm^ 

THIS SUBROUTINE TAKES THE IMAGE POINTS XIMA,YIMA AND 
Oa^VERTS THEM TO I,J SCREEN POINTS 
THIS IS THE AFFIN XPOFM — SEE REF 3 



INRnS = IMAXO 


X image coordinate 


nmY() 


y image coordinate 


5CFRAME 


image coordinate x axis size 


YFERAME 


image coordinate y axis size 


ITOT 


number of elev points 


HIDDENO 


hidden point flag 


OUTEUTS = IA() 


screen x coordinate 


JA() 


screen y coordinate 



****************************************************************** 
IMPLICTT INTEGER (A-Z) 

0CMM:»T /EIEV/ ROW_El^,CX)L_EIiV,TOT_ELEV 

INTEGER IA(TOT_EI£V) , JA(TOT_EIEV) , HIDDEN (TOT_EI£V) ,ITOT 

REAL*8 XIMA,YIMA,A1,A2,B1,B2,C1,C2,DENCM 

REAL*8 IMAX(TOT_EI£V),IMAY(TDT_EIEV)fXFRAME,YFRAME 

DATA I_MAX,J_MAX/512. 0,512.0/ 

calc the affin parameters 



A1=XFRAME/ (I_MAX*1. 0) 

A2=0.0 

B1=0.0 

B2=YFRAME/ ( J_MAX*1 . 0) 

Cl=-XFRAME/2.0 

C2=-YFRAME/2.0 

C 

DO 25 IR=l,rrOT 

IF (HIDDEN(IR) .EQ. 1) GO TO 25 
XIMA=IMAX(IR) 

YIMA=IMAY(IR) 
lA ( IR) = (XIMA-Cl) /A1 
JA ( IR) = ( YIMA-C2 ) /B2 
25 OCmiNUE 

REIURN 
END 
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*************************************************************** 

SUBRCUriNE HIDDEN_SURF(IA,JA,IENEM,IENEN, DEPTH, 

PLJ^, VD_ID , HIDDEN , DEPfflJKEY , WD_ID) 

THIS ROOTINE WILL CHECK EACH PANEL AGAINST ALL 
OTHERS FOR OVERIAP. 

INPUTS => IA( ) 

JA( ) 

lENEN 

niNCM 

HnX)EN() flag for points hidden behind POV 
OOTFUTS = NONE 

*************************************************************** 
IMPLICIT INTEGER (A-Z) 

ocmm:»^ /EIEV/ RCW_EIEV,COL_EIEV,TOT_EIEV,NPLANES 

INTEGER IGREY,NODE_A,NODE_B,NODE_C 
INTEGER IHCN,IENEM,IR,I,J 

INTEGER IPIANES, HIDDEN (TOT_ELEV) ,PL_AR(NPLANES, 5) 

INTEGER lA (TOT_ELEV) , JA (TOT_ELEV) , DEPIH_KEY (NPLANES ) 

INTEGER X1,Y1,X2,Y2,X3,Y3 
REAL*8 DEPTH (TOT_EI£V) 

IJDGICAL SORTED, PIUS, MINUS 

calculate the number of planes 

IPIANES=( (IENEM-1) *2) *(IENEN-1) 

O0UNT=0 

read in the data for each plane 

DO 100 IR=1, IPIANES 
PTA=PL_AR(IR, 1) 

PTB=PL_AR(IR,2) 

PIC=PL_AR(rR, 3) 

IF (HIDDEN(PTA) .EQ. 1 .OR. HIDDEN(PTB) .EQ. 1 .OR. HIDDEN(PTC) 
.EQ. 1) THEN 
GO TO 100 
ELSE 

CX)UNT=OOUNIH-1 
DEPIH_KEY (OOUNT) =IR 
ENDIF 
>0 cx^rnNUE 

sort d^jth values 
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c 

IPLANES==OOUNT 
DO 105 L=1,IPIANES 
IR=DEPIH_KEY(L) 

PLJ^(IR, 5) =MAX(DEPIH(PL_AR(IR, 1 ) ) f 
DEPIH(PLJ^(IR,2) ) ,DEPTH(PL_AR(IR,3) ) ) 

105 OONTINUE 
C 

CALL QUI(3SQRr(PL_AR,IPIANES,DEPlH_KEY) 

C 

CALL UISDC$EE?ASE (WD_ID) 

C 

C now begin selection 
C 

109 DO no K=IPIANES,1,-1 

ir=depih_ke:y(K) 

c 

C draw to virtual di^lay 
C 

NODE_A=PL_AR(IR, 1) 

NODE_B=PL_AR ( IR, 2 ) 

NODE_C=PL_AR (IR, 3 ) 

IGREY=PL_AR(IR,4) 

C 

C limit the gray scale to 250 shades of gray 
C vice 256 shades 
C 

IF (IGREY .EQ. 0) IGREy=l 
IF (IGREY .GT. 250) IGREY=250 
Xl=IA(NODE_A) 

Yl=JA(NODE_A) 

X2=IA(NODE_B) 

Y2<TA(NODE_B) 

X3=IA(NODE_C) 

Y3<TA(NODE_C) 

CALL UIS$SETJWRITIM3_INDEX (VD_ID, 1,1/ IC3^) 
CALL UISDC$PlJ0r(WD_ID,l,Xl,Yl,X3,Y3,X2,Y2) 
C 

no CX^TTINUE 
RETURN 
END 
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********************************************** 
SUBROUTINE CU^CRN 
This routine clears the screen 
Irputs = NONE 
Outputs = NCX® 

********************************************** 
IMPLICIT INTEGER (A-Z) 

CALL SM3$CREATE_PASTEBCARD(PB_ID) 

CALL SMG$ERASE_PASTEBQARD(PB_ID) 

RETURN 

END 



***************************************************** 
SUBROUTINE GC»IV2SEC(DEG,MIN,SEC) 



C ***************************************************** 
IMPLZCTT INTEGER (A-Z) 

IF (DEG .GE. 0) THEN 



SEC=SEC4MIN*6CH-DEG*3600 

EISE 

SEO=DEG*3600-MIN*60-SEC 

ENDIF 

RETURN 

END 



C 

C 

C 

C 

C 

C 

C 



OUTEUTS= SEC 



INKTIS = DEG 



MIN 

SEC 



95 



o o 



************************************************** 



c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



SUBRDOTINE QUICKSORT (ARI^Y, COUNT, KEY) 



INTUT = ARRAYO 
KEY() 
COUNT 



array to sort, sort on 5th field 
key index to main array 
nuniber of elements to sort 



OUTFUT= KEY() new key index to main array 
**************************************************** 



IMPLECTT INTEGER (A-Z) 



OOM^O^/EI^/NN, m, TOT, NPIANES 
INTEGER ST(100,2) 

INTEGER COUNT, KEY (NPIANES) 
INTEGER ARRAY (NPIANES, 5) 



SP=1 

ST(1,1)=1 
ST(1,2)=OOUNT 
DO WHIIE (SP .NE. 0) 
Ir=ST(SP,l) 
R=ST(SP,2) 

SP^P-1 

DO WHILE (R .GT. L) 
LI=L 



RI=R 

INDEX=INT( (L+-R)/2) 

S^^^ARRAY (KEY (INDEX) ,5) 

DO WHIIE (U .IE. RI) 

DO WHILE (ARRAY (KEY (U), 5) .IT. SA) 
LI=LI+1 
END DO 

DO WHILE (ARRAY(KEY(RI) ,5) .GT. SA) 
RI=RI-1 
END DO 

IF (IT .IE. RI) THEN 
TEME^KEY(IT) 

KEY(IT)=KEY(RI) 

KEY(RI)=TE11P 

IT=IT+1 



RI=RI-1 



ENDIF 
END DO 

IF ((R-IT) .GT. (RI-L)) THEN 
IF (L .IT. RI) THEN 
SP=SPfl 
ST(SP,1)=L 
ST(SP,2)=RI 
ENDIF 
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L=LI 

FT.CT F. 

IF (LI .1ST. R) THEN 
SP=SPfl 
ST(SP,1)=LI 
ST(SP,2)=R 
ENDIF 
R=RI 
ENDIF 
END DO 
END DO 
RETURN 
EUD 
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c************************************************* 

c 

SUERXrriNE IMAGE_SAVE(V€)_ID,FII£JIAME, 

C 



c 

c 



IMPLICIT IHrE3GER(A-Z) 



C 



CHARACTER FHE_NAME*(*) ,NAME1*23,NAME2*23 



INTEGER WD_ID,FIAG 
C 

DATA RETLEN/0/ 

DATA RSraTUVO/ 

DATA RHEIGHT/O/ 
DATA BPPIXEI/0/ 



determine the size of the image buffer 

CALL UISDC$READ_IMAGE (WD_ID, 0,0,512,512, IWIDIH, RHEIGHT, BPPIXEL, 
ENOODED1,0) 

RETIEN=F5^DIH*RHEIGHr 

open the external files 

IF (FIAG .EQ. 1) THEN 
C1=INDEX(FII£_NAME, ' . ' ) 

IF (Cl .EQ. 0) THEN 

C2=INDEX(FILE_NAME. ' ') 

NAME1=FII£_NAME ( 1 : C2-1) // ' . SIZ ' 

NAME2=FII£_NAME ( 1 : C2-1) //'.IMS' 

ELSE 

NAME1=FII£_NAME(1:C1-1) // '.SIZ' 

NAME2=FII£_NAME 

ENDIF 

WRITE ( 6 , * ) NAMEl , NAME2 

OPEN (UNIT^11,FIIF>=4^AME1,STAIUS= 'NEW' , 

AOCESS= ' SEQUENTIAL' , PC)FT^= 'UNFORMATTED ' ) 



OPEN (UNIT=10 , FIIE=^^AME2 , STAIUS= 'NEW' , 

F0RM= ' UNFORMATTED ' ) 

WRITE ( 11) RS^IDIH, RHEIGHT, BPPIXEL, RETLEN 

cillocate virtual memory for the buffer 

STATUS=LrB$GET_VM (RETIEN, ENCODED) 

IF (.NOT. STATUS) CALL LrB$STOP(%VAL( STATUS ) ) 
ENDIF 

extract and store private data in buffer 
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c 



CALL UISDC$READ_IMAGE (WD_ID, 0,0,512,512, R^DIH, RHEIGHT, BPPIXEL, 
%VAL (ENCODED) ,PETLEN) 

call subroutine to write the contents of the buffer 

CALL BUFFEK^RTTE (%VAL(ENOODED) ,RETLEN,10) 

RETORN 

END 



SUBROC7ITNE BUFFEf^^RITE (BUE'FER, I£NG?IH, ILN) 
********************************************* 

IMPLICIT INTEGER (A-Z) 

C 

BYTE BUFFER (lENGflH) 

C 

WRITE (ILN) BUFFER 
C 

RETORN 

END 
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